Manuscript Figures
Charlie Bayne
2024-04-05
#Loading in mapping data that will allow us to go from Uniprot ID to gene ID.
# Will use this to fix the labels.
# Next time make sure to have fragpipe give gene level data instead of by uniprot id.
protein_to_gene = readr::read_csv("00_r_input_data/protein_id_to_gene.csv")Figure 1 - Overview of proteomics and metabolomics data
A - Description of the overall study.
This figure was generated in Biorender.
B- Characteristics of patient samples.
### Figure 1b ###
# Reading in the cleaned clinical metadata to start to make plots
clinical_metadata = readr::read_csv("00_r_input_data/clinical_metadata.csv")
p0 = clinical_metadata %>%
ggplot(aes(condition,fill = condition))+
geom_histogram(stat = "count")+
scale_fill_viridis_d()+
ggtitle("Condition")+
ggprism::theme_prism(base_size = 6)
p0
# Defining the columns that we want to show patient charachteristics for
columns_to_look_at = c("gender", "charleson_index", "bacteremia_duration","gender","day_of_blood_draw","sensitive_to_vancomycin","condition","death_during_admission")
# Filtering the data to the relevant columns
sel = clinical_metadata %>%
select(sample_id,columns_to_look_at,pathogen) %>%
filter(!is.na(pathogen))
# Plotting Gender by Pathogen
p1 = sel %>%
ggplot(aes(gender,fill = pathogen))+
geom_histogram(stat = "count")+
scale_fill_viridis_d(end = 0.5)+
ggtitle("Gender") +
ggprism::theme_prism(base_size = 6)
#Plotting Charleson index by Pathogen
p2 = sel %>%
ggplot(aes(charleson_index,fill = pathogen))+
geom_histogram(stat = "count")+
scale_fill_viridis_d(end = 0.5)+
ggtitle("Charleson Index")+
ggprism::theme_prism(base_size = 6)
# Bacteremia Durration
p3 = sel %>%
ggplot(aes(bacteremia_duration,fill = pathogen))+
geom_histogram(stat = "count")+
scale_fill_viridis_d(end = 0.5)+
ggtitle("Bacteremia Duration")+
ggprism::theme_prism(base_size = 6)
# Day of Blood Draw
p4 = sel %>%
ggplot(aes(day_of_blood_draw,fill = pathogen))+
geom_histogram(stat = "count")+
scale_fill_viridis_d(end = 0.5)+
ggtitle("Day of Blood Draw")+
ggprism::theme_prism(base_size = 6)
# Sensitivity to Vancomycin
p5 = sel %>%
ggplot(aes(sensitive_to_vancomycin,fill = pathogen))+
geom_histogram(stat = "count")+
scale_fill_viridis_d(end = 0.5)+
ggtitle("Vancomycin Sensitivity")+
ggprism::theme_prism(base_size = 6)
# Sensitivity to Vancomycin
p6 = sel %>%
ggplot(aes(death_during_admission,fill = pathogen))+
geom_histogram(stat = "count")+
scale_fill_viridis_d(end = 0.5)+
ggtitle("Death During Admission")+
ggprism::theme_prism(base_size = 6)
f1b = ggpubr::ggarrange(p1,p2,p3,p4,p5,p6,common.legend = TRUE,legend = "none",ncol = 6,font.label = list(size = 14))
f1b
ggsave("02_figures/f1_overview_of_study/f1b.pdf",plot = f1b, height = 2.5, width = 8.268,units = "in")C- Overall clustering of proteomoics data.
### Figure 1C ###
#Loading in mapping file
mapping = readxl::read_excel("../proteomics/sample_mapping.xlsx")
# Loading in the clinical metadata
cmd = read_csv("00_r_input_data/clinical_metadata.csv") %>%
select(-condition)
mapping = inner_join(mapping,cmd,by = c("sample_name" = "sample_id"))
#Loading in data
data = read_delim("../proteomics/fragpipe_output_fixed_labels/tmt-report/abundance_gene_MD.tsv")
#Converting to long format
long = data %>%
pivot_longer(6:length(.),names_to = "sample_name")
#Combining data and mapping
all = inner_join(mapping,long,by="sample_name")
#Converting into a matrix
mat = all %>%
dplyr::select(-ProteinID,-NumberPSM,-ReferenceIntensity,-MaxPepProb) %>%
pivot_wider(values_from = value,names_from = Index)
#Extracting only the data
dat = mat %>%
dplyr::select(81:length(.)) %>%
dplyr::select_if(~ !any(is.na(.)))
#Extracting only the metadata
md = mat %>%
dplyr::select(1:81)
#setting row names so we can look at the labels if need be
rownames(dat) <- paste0(md$sample_name)
# Calculating distances
d = dist(dat)
# Showing as a dendrogram
dend = as.dendrogram(hclust(d,method = "ward.D2"))
library(circlize)
library(dendextend)
# adding color by the type of infection
colors_to_use = viridis::viridis(3)[as.numeric(as.factor(md$condition))]
# getting everything in the correct order
colors_to_use = colors_to_use[order.dendrogram(dend)]
# adding mortality
colors_to_use_2 = c("#000000","#FF0000")[as.numeric(as.factor(md$death_during_admission))]
colors_to_use_2 = colors_to_use_2[order.dendrogram(dend)]
# adding color to the dendrogram branches
dend = branches_color(dend,col = colors_to_use) %>%
set("branches_lwd", 8)
labels_colors(dend) = colors_to_use_2
#this is in circular form
pdf("02_figures/f1_overview_of_study/f1c.pdf")
circlize_dendrogram(dend,
labels_track_height = NA,
dend_track_height = 0.5,
labels = TRUE)
dev.off()## png
## 2circlize_dendrogram(dend,
labels_track_height = NA,
dend_track_height = 0.5,
labels = TRUE) 
# stats for dendrogram by k = 2 cluster
circlize_dendrogram(dend %>% color_branches(k = 2),
labels_track_height = NA,
dend_track_height = 0.5,
labels = TRUE) 
D- Overall clustering of meatabolomics data.
### Figure 1 D ###
## Reading in data ##
normalized_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
##Why do we have different numbers of healthy?
test = normalized_data %>%
select(sample_id,condition) %>%
distinct() %>%
group_by(condition) %>%
summarise(n = n())
## Transforming the data ##
norm_wide = normalized_data %>%
select(-row_m_z,-row_retention_time,-value) %>%
tidyr::pivot_wider(names_from = row_id,values_from = norm_value)
norm_dat = norm_wide %>%
dplyr::select(78:length(.))
norm_md = norm_wide %>%
dplyr::select(1:77)
rownames(norm_md) = norm_md$sample_id
## Plotting the overal clustering ##
rownames(norm_dat) = norm_md$sample_id
d = dist(norm_dat)
dend = as.dendrogram(hclust(d, "ward.D2"))
#dend = as.dendrogram(hclust(d))
library(circlize)
library(dendextend)
# let's add some color:
colors_to_use = viridis::viridis(3)[as.numeric(as.factor(norm_md$condition))]
#getting in the same order as dendrogram
colors_to_use = colors_to_use[order.dendrogram(dend)]
# adding mortality
colors_to_use_2 = c("#000000","#FF0000")[as.numeric(as.factor(norm_md$death_during_admission))]
colors_to_use_2 = colors_to_use_2[order.dendrogram(dend)]
# Now we can color the brancehes
dend = branches_color(dend,col = colors_to_use,) %>%
set("branches_lwd", 8)
labels_colors(dend) = colors_to_use_2
pdf("02_figures/f1_overview_of_study/f1d.pdf")
circlize_dendrogram(dend,
labels_track_height = NA,
dend_track_height = 0.5,
labels = TRUE)
dev.off()## png
## 2circlize_dendrogram(dend,
labels_track_height = NA,
dend_track_height = 0.5,
labels = TRUE)
Figure 2 - Host Response Bacteremias Proteomics
A- Volcano plot of healthy vs Infected EcB
Data Processing
# Read MSstats.csv file.
data = read_csv("00_r_input_data/msstats.csv")
#Note that the condition of the bridge channel needs to be labeled as Norm
annotation = read_csv("00_r_input_data/MSstatsTMT_annotation.csv")
mstats = MSstatsTMT::PhilosophertoMSstatsTMTFormat(input = data,annotation)## INFO [2024-04-01 14:44:32] ** Raw data from Philosopher imported successfully.
## INFO [2024-04-01 14:44:34] ** Using provided annotation.
## INFO [2024-04-01 14:44:34] ** Run and Channel labels were standardized to remove symbols such as '.' or '%'.
## INFO [2024-04-01 14:44:34] ** The following options are used:
## - Features will be defined by the columns: PeptideSequence, PrecursorCharge
## - Shared peptides will be removed.
## - Proteins with single feature will not be removed.
## - Features with less than 3 measurements within each run will be removed.
## INFO [2024-04-01 14:44:34] ** Rows with values not greater than 0.6 in Purity are removed
## WARN [2024-04-01 14:44:34] ** PeptideProphetProbability not found in input columns.
## INFO [2024-04-01 14:44:34] ** Sequences containing Oxidation are removed.
## INFO [2024-04-01 14:44:35] ** Features with all missing measurements across channels within each run are removed.
## INFO [2024-04-01 14:44:36] ** Shared peptides are removed.
## INFO [2024-04-01 14:44:36] ** Features with one or two measurements across channels within each run are removed.
## INFO [2024-04-01 14:45:52] ** PSMs have been aggregated to peptide ions.
## INFO [2024-04-01 14:45:54] ** Run annotation merged with quantification data.
## INFO [2024-04-01 14:45:56] ** For peptides overlapped between fractions of plex1_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:45:56] ** For peptides overlapped between fractions of plex2_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:45:56] ** For peptides overlapped between fractions of plex3_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:45:56] ** For peptides overlapped between fractions of plex4_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:45:56] ** For peptides overlapped between fractions of plex5_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:45:56] ** For peptides overlapped between fractions of plex6_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:45:56] ** For peptides overlapped between fractions of plex7_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:45:57] ** Fractions belonging to same mixture have been combined.
## INFO [2024-04-01 14:45:57] ** Features with one or two measurements across channels within each run are removed.
## INFO [2024-04-01 14:45:57] ** Fractionation handled.
## INFO [2024-04-01 14:45:57] ** Updated quantification data to make balanced design. Missing values are marked by NA
## INFO [2024-04-01 14:45:57] ** Finished preprocessing. The dataset is ready to be processed by the proteinSummarization function.stats = MSstatsTMT::PhilosophertoMSstatsTMTFormat(input = data,annotation)## INFO [2024-04-01 14:45:57] ** Raw data from Philosopher imported successfully.
## INFO [2024-04-01 14:45:59] ** Using provided annotation.
## INFO [2024-04-01 14:45:59] ** Run and Channel labels were standardized to remove symbols such as '.' or '%'.
## INFO [2024-04-01 14:45:59] ** The following options are used:
## - Features will be defined by the columns: PeptideSequence, PrecursorCharge
## - Shared peptides will be removed.
## - Proteins with single feature will not be removed.
## - Features with less than 3 measurements within each run will be removed.
## INFO [2024-04-01 14:45:59] ** Rows with values not greater than 0.6 in Purity are removed
## WARN [2024-04-01 14:45:59] ** PeptideProphetProbability not found in input columns.
## INFO [2024-04-01 14:45:59] ** Sequences containing Oxidation are removed.
## INFO [2024-04-01 14:46:00] ** Features with all missing measurements across channels within each run are removed.
## INFO [2024-04-01 14:46:01] ** Shared peptides are removed.
## INFO [2024-04-01 14:46:01] ** Features with one or two measurements across channels within each run are removed.
## INFO [2024-04-01 14:47:17] ** PSMs have been aggregated to peptide ions.
## INFO [2024-04-01 14:47:18] ** Run annotation merged with quantification data.
## INFO [2024-04-01 14:47:20] ** For peptides overlapped between fractions of plex1_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:47:20] ** For peptides overlapped between fractions of plex2_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:47:21] ** For peptides overlapped between fractions of plex3_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:47:21] ** For peptides overlapped between fractions of plex4_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:47:21] ** For peptides overlapped between fractions of plex5_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:47:21] ** For peptides overlapped between fractions of plex6_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:47:21] ** For peptides overlapped between fractions of plex7_1 use the fraction with maximal average abundance.
## INFO [2024-04-01 14:47:21] ** Fractions belonging to same mixture have been combined.
## INFO [2024-04-01 14:47:21] ** Features with one or two measurements across channels within each run are removed.
## INFO [2024-04-01 14:47:21] ** Fractionation handled.
## INFO [2024-04-01 14:47:22] ** Updated quantification data to make balanced design. Missing values are marked by NA
## INFO [2024-04-01 14:47:22] ** Finished preprocessing. The dataset is ready to be processed by the proteinSummarization function.# use MSstats for protein summarization
quant.msstats = MSstatsTMT::proteinSummarization(stats,
method="msstats",
global_norm=TRUE,
reference_norm=TRUE,
remove_norm_channel = TRUE,
remove_empty_channel = TRUE)## INFO [2024-04-01 14:47:22] ** MSstatsTMT - proteinSummarization function
## INFO [2024-04-01 14:47:22] Summarizing for Run : plex4_1 ( 1 of 7 )
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
## INFO [2024-04-01 14:47:37] Summarizing for Run : plex2_1 ( 2 of 7 )
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
## INFO [2024-04-01 14:47:56] Summarizing for Run : plex6_1 ( 3 of 7 )
## | | | 0% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%
## INFO [2024-04-01 14:48:07] Summarizing for Run : plex7_1 ( 4 of 7 )
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
## INFO [2024-04-01 14:48:17] Summarizing for Run : plex3_1 ( 5 of 7 )
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
## INFO [2024-04-01 14:48:33] Summarizing for Run : plex5_1 ( 6 of 7 )
## | | | 0% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%
## INFO [2024-04-01 14:48:45] Summarizing for Run : plex1_1 ( 7 of 7 )
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
## INFO [2024-04-01 14:49:05] ** Protein-level summarization done by MSstats.# Taking the pl_data and formatting appropriately
pl_data = quant.msstats$ProteinLevelData %>%
tidyr::separate(Protein, c("orientation","uniprotid","locus"),sep ="\\|") %>%
# Converting to gene name
inner_join(protein_to_gene,by = c("uniprotid" = "ProteinID")) %>%
mutate(sample_id = paste0(Mixture,".",Channel))
write_csv(pl_data,"01_r_processed_data/quant_msstats_protein_level_data.csv")
write_csv(quant.msstats$FeatureLevelData,"01_r_processed_data/quant_msstats_feature_level_data.csv")
### Preforming the statistics ###
test.pairwise = MSstatsTMT::groupComparisonTMT(quant.msstats, moderated = TRUE)## INFO [2024-04-01 14:49:06] Design: 7 mixtures.
## INFO [2024-04-01 14:49:06] Design: 1 MS run per mixture.
## INFO [2024-04-01 14:49:06] Design: group comparison design (Different conditions contains different biological subjects).
## INFO [2024-04-01 14:49:06] Model fitting for 598 proteins.
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
## INFO [2024-04-01 14:49:16] Testing for 598 proteins:
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%## All groups ##
res = test.pairwise$ComparisonResult %>%
separate(Protein, c("orientation","uniprotid","locus"),sep ="\\|")
write_csv(res,"01_r_processed_data/groupComparison_results.csv")
## Faecalis v Healthy
faecalis_v_healthy = res %>%
filter(Label == "Faecalis vs Healthy")
write_csv(faecalis_v_healthy,"01_r_processed_data/faecalis_v_healthy.csv")
## Faecium v Healthy
healthy_v_faecium = res %>%
filter(Label == "Healthy vs Faecium")
write_csv(healthy_v_faecium,"01_r_processed_data/healthy_v_faecium.csv")
## Faecalis v Faecium ##
faecalis_v_faecium = res %>%
filter(Label == "Faecalis vs Faecium")
write_csv(faecalis_v_faecium,"01_r_processed_data/faecalis_v_faecium.csv")How many are differentially abundant?
infected_vs_healthy = readr::read_csv("01_r_processed_data/Infected_vs_Healthy.csv") %>%
left_join(protein_to_gene,by = c("uniprotid"= "ProteinID"))
infected_vs_healthy %>%
mutate(sign = case_when(log2FC < 0 ~ "negative",
TRUE ~ "positive")) %>%
filter(adj.pvalue <= 0.05) %>%
group_by(sign) %>%
summarise(n = n())infected_vs_healthy = readr::read_csv("01_r_processed_data/Infected_vs_Healthy.csv")
faecalis_v_healthy = readr::read_csv("01_r_processed_data/faecalis_v_healthy.csv")
healthy_v_faecium = readr::read_csv("01_r_processed_data/healthy_v_faecium.csv")
both_sig = infected_vs_healthy %>%
filter(adj.pvalue <= 0.05) %>%
pull(uniprotid)
# Defining significant proteins in faecalis
faecalis_sig = faecalis_v_healthy %>%
filter(adj.pvalue <= 0.05) %>%
pull(uniprotid)
# Defining significant proteins in faecium
faecium_sig = healthy_v_faecium %>%
filter(adj.pvalue <= 0.05) %>%
pull(uniprotid) B - Heat map by GO Terms
Processing/ Generation of related supplementary figures
EcB
Up.
Network analysis of proteins that are shared between faecalis and faecium relative to healthy when considered alone and increased in infected samples
##Loading all the go stuff
#Getting all the go annotations
go = read_csv("00_r_input_data/go_annotations.csv") %>%
dplyr::select(1,length(.)) %>%
tidyr::separate_rows(Gene.Ontology.IDs, sep = ';') %>%
dplyr::mutate(Gene.Ontology.IDs = gsub(" ","",Gene.Ontology.IDs)) %>%
# Going from Protein ID to gene ID
left_join(protein_to_gene,by = c("Entry" = "ProteinID")) %>%
select(-Entry) %>%
select(Entry = Gene,Gene.Ontology.IDs)
#usng enricher
TERM2GENE = go %>%
dplyr::select(TERM = Gene.Ontology.IDs,GENE = Entry) %>%
as.data.frame()
TERM2NAME = AnnotationDbi::select(GO.db::GO.db,keys = unique(TERM2GENE$TERM),columns = c("TERM")) %>%
as.data.frame()
# Extracting the uniprotids up in infected and shared
infected_vs_healthy = readr::read_csv("01_r_processed_data/Infected_vs_Healthy.csv") %>%
left_join(protein_to_gene,by = c("uniprotid"= "ProteinID"))
shared_up_in_infected = infected_vs_healthy %>%
filter(log2FC > 0) %>%
filter(adj.pvalue <= 0.05) %>%
pull(Gene)
# Perfomring GO enrichment
enrich = clusterProfiler::enricher(shared_up_in_infected,
universe = go$Entry,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME,
maxGSSize = 1000)
DT::datatable(enrich@result)EcB_enrich_up = clusterProfiler::enricher(shared_up_in_infected,
universe = go$Entry,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME)
up_ECB = clusterProfiler::cnetplot(enrich,
categorySize="pvalue",
circular = FALSE,
categorySize="pvalue",
cex_label_category = 0.9,
showCategory = 100,
cex_label_gene = 0.5
)+
theme(legend.position = "none")
up_ECB
ggsave("03_supplemental/s2_networks_compared_to_healthy/up_ECB.pdf",up_ECB,height=10,width=10)Down
##Loading all the go stuff
#Getting all the go annotations
go = read_csv("00_r_input_data/go_annotations.csv") %>%
dplyr::select(1,length(.)) %>%
tidyr::separate_rows(Gene.Ontology.IDs, sep = ';') %>%
dplyr::mutate(Gene.Ontology.IDs = gsub(" ","",Gene.Ontology.IDs)) %>%
# Going from Protein ID to gene ID
left_join(protein_to_gene,by = c("Entry" = "ProteinID")) %>%
select(-Entry) %>%
select(Entry = Gene,Gene.Ontology.IDs)
#usng enricher
TERM2GENE = go %>%
dplyr::select(TERM = Gene.Ontology.IDs,GENE = Entry) %>%
as.data.frame()
TERM2NAME = AnnotationDbi::select(GO.db::GO.db,keys = unique(TERM2GENE$TERM),columns = c("TERM")) %>%
as.data.frame()
# Extracting the uniprotids up in infected and shared
shared_down_in_infected = infected_vs_healthy %>%
filter(log2FC < 0) %>%
filter(adj.pvalue <= 0.05) %>%
pull(Gene)
# Perfomring GO enrichment
enrich = clusterProfiler::enricher(shared_down_in_infected,
universe = go$Entry,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME,
maxGSSize = 1000)
DT::datatable(enrich@result)down_EcB = clusterProfiler::cnetplot(enrich,
categorySize="pvalue",
cex_label_category = 1,
showCategory = 100,
cex_label_gene = 0.5)+
theme(legend.position = "none")
down_EcB
ggsave("03_supplemental/s2_networks_compared_to_healthy/down_EcB.pdf",down_EcB,height =10, width = 10)Faecalis
Up
#usng enricher
TERM2GENE = go %>%
dplyr::select(TERM = Gene.Ontology.IDs,GENE = Entry) %>%
as.data.frame()
TERM2NAME = AnnotationDbi::select(GO.db::GO.db,keys = unique(TERM2GENE$TERM),columns = c("TERM")) %>%
as.data.frame()
faecalis_v_healthy = readr::read_csv("01_r_processed_data/faecalis_v_healthy.csv") %>%
left_join(protein_to_gene,by = c("uniprotid"= "ProteinID"))
# Extracting the uniprotids up in infected and shared
up_in_faecalis = faecalis_v_healthy %>%
filter(log2FC >= 0) %>%
filter(adj.pvalue <= 0.05) %>%
pull(Gene)
# Perfomring GO enrichment
enrich = clusterProfiler::enricher(up_in_faecalis,
universe = go$Entry,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME,
maxGSSize = 1000)
DT::datatable(enrich@result)up_faecalis = clusterProfiler::cnetplot(enrich,
categorySize="pvalue",
circular = FALSE,
categorySize="pvalue",
cex_label_category = 1,
showCategory = 100,
cex_label_gene = 0.5)+
theme(legend.position = "none")
up_faecalis
ggsave("03_supplemental/s2_networks_compared_to_healthy/up_faecalis.pdf",up_faecalis,height=10,width=10)Down
#usng enricher
TERM2GENE = go %>%
dplyr::select(TERM = Gene.Ontology.IDs,GENE = Entry) %>%
as.data.frame()
TERM2NAME = AnnotationDbi::select(GO.db::GO.db,keys = unique(TERM2GENE$TERM),columns = c("TERM")) %>%
as.data.frame()
# Extracting the uniprotids up in infected and shared
down_in_faecalis = faecalis_v_healthy %>%
filter(log2FC <= 0) %>%
filter(adj.pvalue <= 0.05) %>%
pull(Gene)
# Perfomring GO enrichment
enrich = clusterProfiler::enricher(down_in_faecalis,
universe = go$Entry,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME,
maxGSSize = 1000)
DT::datatable(enrich@result)down_faecalis = clusterProfiler::cnetplot(enrich,
categorySize="pvalue",
circular = FALSE,
categorySize="pvalue",
cex_label_category = 1,
showCategory = 100,
cex_label_gene = 0.5)+
theme(legend.position = "none")
down_faecalis
ggsave("03_supplemental/s2_networks_compared_to_healthy/down_faecalis.pdf",down_faecalis,height=10,width=10)Faecium
Up
healthy_v_faecium = readr::read_csv("01_r_processed_data/healthy_v_faecium.csv") %>%
left_join(protein_to_gene,by = c("uniprotid"= "ProteinID"))
#usng enricher
TERM2GENE = go %>%
dplyr::select(TERM = Gene.Ontology.IDs,GENE = Entry) %>%
as.data.frame()
TERM2NAME = AnnotationDbi::select(GO.db::GO.db,keys = unique(TERM2GENE$TERM),columns = c("TERM")) %>%
as.data.frame()
# Extracting the uniprotids up in infected and shared
# Extracting the uniprotids up in infected and shared
up_in_faecium = healthy_v_faecium %>%
filter(log2FC <= 0) %>%
filter(adj.pvalue <= 0.05) %>%
pull(Gene)
# Perfomring GO enrichment
enrich = clusterProfiler::enricher(up_in_faecium,
universe = go$Entry,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME,
maxGSSize = 1000)
DT::datatable(enrich@result)up_faecium = clusterProfiler::cnetplot(enrich,
categorySize="pvalue",
circular = FALSE,
categorySize="pvalue",
cex_label_category = 1,
showCategory = 100,
cex_label_gene = 0.5)+
theme(legend.position = "none")
up_faecium
ggsave("03_supplemental/s2_networks_compared_to_healthy/up_faecium.pdf",up_faecium,height=10,width=10)Down
#usng enricher
TERM2GENE = go %>%
dplyr::select(TERM = Gene.Ontology.IDs,GENE = Entry) %>%
as.data.frame()
TERM2NAME = AnnotationDbi::select(GO.db::GO.db,keys = unique(TERM2GENE$TERM),columns = c("TERM")) %>%
as.data.frame()
# Extracting the uniprotids up in infected and shared
down_in_faecium = healthy_v_faecium %>%
filter(log2FC >= 0) %>%
filter(adj.pvalue <= 0.05) %>%
pull(Gene)
# Perfomring GO enrichment
enrich = clusterProfiler::enricher(down_in_faecium,
universe = go$Entry,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME,
maxGSSize = 1000)
DT::datatable(enrich@result)down_faecium = clusterProfiler::cnetplot(enrich,
categorySize="pvalue",
circular = FALSE,
categorySize="pvalue",
cex_label_category = 1,
showCategory = 100,
cex_label_gene = 0.5)+
theme(legend.position = "none")
down_faecium
ggsave("03_supplemental/s2_networks_compared_to_healthy/down_faecium.pdf",down_faecium,height=10,width=10)Staph Aureus
#Loading in data, getting GO annotations
staph_data = readxl::read_excel("00_r_input_data/staph_proteomics.xlsx",skip = 1)
ids = staph_data %>%
pull(Protein)
# Getting GO annotations by hitting Uniprots API
anot_go = GLabR::annotate_uniprot_single(ids,columns = "accession,id,gene_primary,go,go_c,go_p,go_f,go_id")
write_csv(anot_go,"01_r_processed_data/staph_go_annotations.csv")Up
#Loading in data,
staph_data = readxl::read_excel("00_r_input_data/staph_proteomics.xlsx",skip = 1)
anot_go = read_csv("01_r_processed_data/staph_go_annotations.csv") %>%
select(Entry = Gene.Names..primary.,Gene.Ontology.IDs) %>%
dplyr::select(1,length(.)) %>%
tidyr::separate_rows(Gene.Ontology.IDs, sep = ';') %>%
dplyr::mutate(Gene.Ontology.IDs = gsub(" ","",Gene.Ontology.IDs))
# Going from Protein ID to gene ID
# left_join(protein_to_gene,by = c("Entry" = "ProteinID")) %>%
sig_up = staph_data %>%
filter(Infection_pvalue <= 0.05) %>%
filter(Infection_FC > 1) %>%
pull(Gene)
sig_down = staph_data %>%
filter(Infection_pvalue <= 0.05) %>%
filter(Infection_FC < 1) %>%
pull(Gene)
background = staph_data %>%
pull(Gene)
#usng enricher
TERM2GENE = anot_go %>%
dplyr::select(TERM = Gene.Ontology.IDs,GENE = Entry) %>%
as.data.frame()
test = TERM2GENE %>%
filter(TERM == "GO:0070062")
TERM2NAME = AnnotationDbi::select(GO.db::GO.db,keys = unique(TERM2GENE$TERM),columns = c("TERM")) %>%
as.data.frame()
# Perfomring GO enrichment
enrich = clusterProfiler::enricher(sig_up,
universe = background,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME,
maxGSSize = 1000)
DT::datatable(enrich@result)up_staph = clusterProfiler::cnetplot(enrich,
categorySize="pvalue",
cex_label_category = 1,
showCategory = 100,
cex_label_gene = 0.5)+
theme(legend.position = "none")
up_staph
ggsave("03_supplemental/s2_networks_compared_to_healthy/up_staph.pdf",up_staph,height=10,width=10)Down
# Perfomring GO enrichment
enrich = clusterProfiler::enricher(sig_down,
universe = background,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME,
maxGSSize = 1000)
DT::datatable(enrich@result)# Perfomring GO enrichment
staph_down_enrich = clusterProfiler::enricher(sig_down,
universe = background,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME)
down_staph = clusterProfiler::cnetplot(enrich,
categorySize="pvalue",
cex_label_category = 1,
showCategory = 100,
cex_label_gene = 0.5)+
theme(legend.position = "none")
down_staph
ggsave("03_supplemental/s2_networks_compared_to_healthy/down_staph.pdf",down_staph,height=10,width=10)Plotting Heatmap
## Enterococcal Bacteremia
up_in_EcB_data = up_ECB$data[1:32,] %>%
mutate(type = "EcB")
down_in_EcB_data = down_EcB$data[1:16,] %>%
mutate(size = size * -1) %>%
mutate(type = "EcB")
## E. faecalis
up_in_faecalis_data = up_faecalis$data[1:24,] %>%
mutate(type = "E. faecalis")
down_in_faecalis_data = down_faecalis$data[1:24,] %>%
mutate(type = "E. faecalis") %>%
mutate(size = size * -1)
## E. faecium
up_in_faecium_data = up_faecium$data[1:9,] %>%
mutate(type = "E. faecium")
down_in_faecium_data = down_faecium$data[1:8,] %>%
mutate(type = "E. faecium") %>%
mutate(size = size * -1)
## Staph Aureus
up_staph_data = up_staph$data[1:8,] %>%
mutate(type = "S. aureus")
down_staph_data = down_staph$data[1:19,] %>%
mutate(size = size * -1) %>%
mutate(type = "S. aureus")
# Combining all of the data
all = bind_rows(up_in_EcB_data,
down_in_EcB_data,
up_in_faecalis_data,
down_in_faecalis_data,
up_in_faecium_data,
down_in_faecium_data,
up_staph_data,
down_staph_data)
# plotting the heatmap
f2b = all %>%
ggplot(aes(reorder(name,size),type,fill = size))+
geom_tile()+
coord_flip()+
scale_fill_gradient2(
name = waiver(),
high = "red",
low = "green",
mid = "white",
midpoint = 0,
space = "Lab",
na.value = "grey50",
guide = "colourbar",
aesthetics = "fill"
) +
xlab("GO term")+
ylab("Bacteremia Type")+
ggprism::theme_prism(base_size = 10)+
theme(legend.position = "top")
f2b
ggsave("02_figures/f2_infection_v_healthy_proteomics/f2b.pdf",height = 10, width = 7)? Platelet Alpha Graunle
source("functions.R")
# Extracting the proteins with the term platelet alphaa granule lumen from staph enrichment results
staph_enrich_platelet = staph_down_enrich@result %>%
filter(Description == "platelet alpha granule lumen")
down_in_staph = strsplit(staph_enrich_platelet$geneID,"/")[[1]]
# Ex tracting the proteins with the term alph granule lumen from EcB enrichment results
EcB_enrich_platelet = EcB_enrich_up@result %>%
filter(Description == "platelet alpha granule lumen")
up_in_EcB = strsplit(EcB_enrich_platelet$geneID,"/")[[1]]
# These are the shared proteins at the intersection.
shared = intersect(down_in_staph,up_in_EcB)
shared## [1] "PPBP" "MMRN1" "PF4" "SPARC" "THBS1"# plotting differences to healthy - staph
staph = staph_data %>%
filter(Gene %in% down_in_staph) %>%
select(Gene,NN1:length(.)) %>%
pivot_longer(2:length(.),names_to = "sample_id",values_to = "Abundance") %>%
mutate(Condition = case_when(grepl("HS.*|HM.*",sample_id) ~ "Infected",
TRUE ~ "Healthy"),
Abundance = as.numeric(Abundance))
make_violin_plot_proteomics(staph,shared,ncol = length(shared),pname = "p")
## TableGrob (1 x 5) "arrange": 5 grobs
## z cells name grob
## PPBP 1 (1-1,1-1) arrange gtable[layout]
## MMRN1 2 (1-1,2-2) arrange gtable[layout]
## PF4 3 (1-1,3-3) arrange gtable[layout]
## SPARC 4 (1-1,4-4) arrange gtable[layout]
## THBS1 5 (1-1,5-5) arrange gtable[layout]# plotting differnces to healthy EcB
protein_sum = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
f = protein_sum %>%
filter(Gene %in% shared) %>%
mutate(Condition = case_when(Condition == "Healthy" ~ "Healthy",
TRUE ~ "Infected"))
make_violin_plot_proteomics(f,shared,ncol = length(shared),pname = "p.signif")
## TableGrob (1 x 5) "arrange": 5 grobs
## z cells name grob
## PPBP 1 (1-1,1-1) arrange gtable[layout]
## MMRN1 2 (1-1,2-2) arrange gtable[layout]
## PF4 3 (1-1,3-3) arrange gtable[layout]
## SPARC 4 (1-1,4-4) arrange gtable[layout]
## THBS1 5 (1-1,5-5) arrange gtable[layout]C - Venn Diagram
Up in infected
infected_vs_healthy = readr::read_csv("01_r_processed_data/Infected_vs_Healthy.csv")
faecalis_v_healthy = readr::read_csv("01_r_processed_data/faecalis_v_healthy.csv")
healthy_v_faecium = readr::read_csv("01_r_processed_data/healthy_v_faecium.csv")
# Defining significant proteins in faecalis
faecalis_sig = faecalis_v_healthy %>%
filter(adj.pvalue <= 0.05) %>%
filter(log2FC<0) %>%
pull(uniprotid)
# Defining significant proteins in faecium
faecium_sig = healthy_v_faecium %>%
filter(adj.pvalue <= 0.05) %>%
filter(log2FC > 0) %>%
pull(uniprotid)
# found to be sig in both faecalis and faecium
enterococcus_sig = infected_vs_healthy %>%
filter(adj.pvalue <= 0.05) %>%
filter(log2FC < 0)%>%
pull(uniprotid)
# Loading staph_data
staph_sig = readxl::read_excel("00_r_input_data/staph_proteomics.xlsx",skip = 1) %>%
filter(Infection_pvalue <= 0.05,
Infection_FC < 1) %>%
pull(Protein)
library(ggVennDiagram)
up = ggVennDiagram(list(EcB = enterococcus_sig, faecalis = faecalis_sig, faecium = faecium_sig, S.aureus = staph_sig))+
theme(legend.position = "none")+
scale_fill_distiller(palette = "Reds", direction = 1)
up
ggsave("02_figures/f2_infection_v_healthy_proteomics/f2d_a.pdf", up, height = 5, width = 6)
ggVennDiagram(list(EcB = enterococcus_sig, S.aureus = staph_sig))
Up in Healthy
infected_vs_healthy = readr::read_csv("01_r_processed_data/Infected_vs_Healthy.csv")
faecalis_v_healthy = readr::read_csv("01_r_processed_data/faecalis_v_healthy.csv")
healthy_v_faecium = readr::read_csv("01_r_processed_data/healthy_v_faecium.csv")
# Defining significant proteins in faecalis
faecalis_sig = faecalis_v_healthy %>%
filter(adj.pvalue <= 0.05) %>%
filter(log2FC>0) %>%
pull(uniprotid)
# Defining significant proteins in faecium
faecium_sig = healthy_v_faecium %>%
filter(adj.pvalue <= 0.05) %>%
filter(log2FC < 0) %>%
pull(uniprotid)
# found to be sig in both faecalis and faecium
enterococcus_sig = infected_vs_healthy %>%
filter(adj.pvalue <= 0.05) %>%
filter(log2FC > 0)%>%
pull(uniprotid)
# Loading staph_data
staph_sig = readxl::read_excel("00_r_input_data/staph_proteomics.xlsx",skip = 1) %>%
filter(Infection_pvalue <= 0.05,
Infection_FC > 1) %>%
pull(Protein)
library(ggVennDiagram)
down = ggVennDiagram(list(EcB = enterococcus_sig, faecalis = faecalis_sig, faecium = faecium_sig, S.aureus = staph_sig))+
theme(legend.position = "none")+
scale_fill_distiller(palette = "Greens", direction = 1)
down
ggsave("02_figures/f2_infection_v_healthy_proteomics/f2d_b.pdf",down,height = 5, width = 6)
ggVennDiagram(list(EcB = enterococcus_sig, S.aureus = staph_sig))
What proteins have differing responses in staph vs EcB?
# found to be sig in both faecalis and faecium
enterococcus_sig_up = infected_vs_healthy %>%
filter(adj.pvalue <= 0.05) %>%
filter(log2FC > 0)%>%
pull(uniprotid)
# Loading staph_data
staph_sig_down = readxl::read_excel("00_r_input_data/staph_proteomics.xlsx",skip = 1) %>%
filter(Infection_pvalue <= 0.05,
Infection_FC < 1) %>%
pull(Protein)
library(ggVennDiagram)
ggVennDiagram(list(enterococcus_sig_up, staph_sig_down))
GLabR::annotate_uniprot(intersect(enterococcus_sig_up,staph_sig_down),"accession,id,gene_primary,protein_name")## | | | 0% | |==================================================| 100%# found to be sig in both faecalis and faecium
enterococcus_sig_down = infected_vs_healthy %>%
filter(adj.pvalue <= 0.05) %>%
filter(log2FC < 0)%>%
pull(uniprotid)
# Loading staph_data
staph_sig_up = readxl::read_excel("00_r_input_data/staph_proteomics.xlsx",skip = 1) %>%
filter(Infection_pvalue <= 0.05,
Infection_FC > 1) %>%
pull(Protein)
library(ggVennDiagram)
ggVennDiagram(list(enterococcus_sig_down, staph_sig_up))
GLabR::annotate_uniprot(intersect(enterococcus_sig_down,staph_sig_up),"accession,id,gene_primary,protein_name")## | | | 0% | |==================================================| 100%C- Evaluation of Proteomic Biomarkers
Generating top 10 biomarkers
#reading in the data
pl = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
## Getting into a format ot use with EFS ##
clean = pl %>%
dplyr::group_by(Mixture,Gene,Channel,Condition) %>%
dplyr::summarise(sum = sum(Abundance)) %>%
ungroup()
infected_healthy_l = clean %>%
mutate(classlabels = case_when(Condition %in% c("Faecalis","Faecium") ~ 1,
TRUE ~ 0)) %>%
mutate(classlabels = as.logical(classlabels)) %>%
# - messes up EFS for some reason, very annoying.
mutate(Gene = gsub("-","_",Gene))
infected_healthy_w = infected_healthy_l %>%
tidyr::pivot_wider(names_from = Gene,
values_from = sum) %>%
mutate(sampleid = paste0(Mixture,".",Channel)) %>%
dplyr::select(-Mixture,-Channel,-Condition,-sampleid) %>%
dplyr::select_if(~ !any(is.na(.))) %>%
#needs to be a dataframe or EFS freaks out.
as.data.frame()
set.seed(2)
## Running EFS ##
ensemble_result = EFS::ensemble_fs(infected_healthy_w,1,cor_threshold = 0)## [1] "default value for NA_threshold = 0.2"
## [1] "default value for runs = 100"
## [1] "default value for selection is c(TRUE, TRUE, TRUE,TRUE, TRUE, TRUE, FALSE, FALSE)"
## [1] "Start Median"
## [1] "Start Pearson"
## [1] "Start Spearman"
## [1] "Start LogReg"
## [1] "Start RF"
## [1] 1
## Time difference of 0.2842882 secs
## [1] 2
## Time difference of 0.2777565 secs
## [1] 3
## Time difference of 0.2796893 secs
## [1] 4
## Time difference of 0.2805221 secs
## [1] 5
## Time difference of 0.3027546 secs
## [1] 6
## Time difference of 0.2823405 secs
## [1] 7
## Time difference of 0.2778196 secs
## [1] 8
## Time difference of 0.2801352 secs
## [1] 9
## Time difference of 0.2791164 secs
## [1] 10
## Time difference of 0.372947 secs
## [1] 11
## Time difference of 0.2784245 secs
## [1] 12
## Time difference of 0.2800491 secs
## [1] 13
## Time difference of 0.2758179 secs
## [1] 14
## Time difference of 0.2792733 secs
## [1] 15
## Time difference of 0.2762392 secs
## [1] 16
## Time difference of 0.2804942 secs
## [1] 17
## Time difference of 0.2757726 secs
## [1] 18
## Time difference of 0.2796123 secs
## [1] 19
## Time difference of 0.2788227 secs
## [1] 20
## Time difference of 0.2809844 secs
## [1] 21
## Time difference of 0.2764399 secs
## [1] 22
## Time difference of 0.2813222 secs
## [1] 23
## Time difference of 0.2765057 secs
## [1] 24
## Time difference of 0.2802126 secs
## [1] 25
## Time difference of 0.2790532 secs
## [1] 26
## Time difference of 0.2787097 secs
## [1] 27
## Time difference of 0.2783546 secs
## [1] 28
## Time difference of 0.2798426 secs
## [1] 29
## Time difference of 0.2755103 secs
## [1] 30
## Time difference of 0.2801824 secs
## [1] 31
## Time difference of 0.2795572 secs
## [1] 32
## Time difference of 0.2781978 secs
## [1] 33
## Time difference of 0.277921 secs
## [1] 34
## Time difference of 0.2785676 secs
## [1] 35
## Time difference of 0.2768929 secs
## [1] 36
## Time difference of 0.280632 secs
## [1] 37
## Time difference of 0.2771516 secs
## [1] 38
## Time difference of 0.2786884 secs
## [1] 39
## Time difference of 0.2792988 secs
## [1] 40
## Time difference of 0.2756121 secs
## [1] 41
## Time difference of 0.2809896 secs
## [1] 42
## Time difference of 0.281424 secs
## [1] 43
## Time difference of 0.2811427 secs
## [1] 44
## Time difference of 0.278234 secs
## [1] 45
## Time difference of 0.2796485 secs
## [1] 46
## Time difference of 0.2767174 secs
## [1] 47
## Time difference of 0.2792492 secs
## [1] 48
## Time difference of 0.2793458 secs
## [1] 49
## Time difference of 0.2811592 secs
## [1] 50
## Time difference of 0.2765591 secs
## [1] 51
## Time difference of 0.2807231 secs
## [1] 52
## Time difference of 0.276762 secs
## [1] 53
## Time difference of 0.2800379 secs
## [1] 54
## Time difference of 0.2790682 secs
## [1] 55
## Time difference of 0.3704154 secs
## [1] 56
## Time difference of 0.2788606 secs
## [1] 57
## Time difference of 0.2781799 secs
## [1] 58
## Time difference of 0.2797458 secs
## [1] 59
## Time difference of 0.2799327 secs
## [1] 60
## Time difference of 0.2768764 secs
## [1] 61
## Time difference of 0.2773247 secs
## [1] 62
## Time difference of 0.2804031 secs
## [1] 63
## Time difference of 0.2752573 secs
## [1] 64
## Time difference of 0.2749345 secs
## [1] 65
## Time difference of 0.2795067 secs
## [1] 66
## Time difference of 0.2757018 secs
## [1] 67
## Time difference of 0.2781987 secs
## [1] 68
## Time difference of 0.2767062 secs
## [1] 69
## Time difference of 0.2767646 secs
## [1] 70
## Time difference of 0.2803757 secs
## [1] 71
## Time difference of 0.2818537 secs
## [1] 72
## Time difference of 0.2773337 secs
## [1] 73
## Time difference of 0.2751365 secs
## [1] 74
## Time difference of 0.2802444 secs
## [1] 75
## Time difference of 0.277076 secs
## [1] 76
## Time difference of 0.277879 secs
## [1] 77
## Time difference of 0.2793305 secs
## [1] 78
## Time difference of 0.2770395 secs
## [1] 79
## Time difference of 0.281132 secs
## [1] 80
## Time difference of 0.2782218 secs
## [1] 81
## Time difference of 0.277498 secs
## [1] 82
## Time difference of 0.2797217 secs
## [1] 83
## Time difference of 0.2754848 secs
## [1] 84
## Time difference of 0.2786725 secs
## [1] 85
## Time difference of 0.2784717 secs
## [1] 86
## Time difference of 0.2760696 secs
## [1] 87
## Time difference of 0.2761865 secs
## [1] 88
## Time difference of 0.2801964 secs
## [1] 89
## Time difference of 0.2770946 secs
## [1] 90
## Time difference of 0.2782071 secs
## [1] 91
## Time difference of 0.2772059 secs
## [1] 92
## Time difference of 0.2784343 secs
## [1] 93
## Time difference of 0.278909 secs
## [1] 94
## Time difference of 0.2811127 secs
## [1] 95
## Time difference of 0.2760377 secs
## [1] 96
## Time difference of 0.2776377 secs
## [1] 97
## Time difference of 0.2788477 secs
## [1] 98
## Time difference of 0.2805719 secs
## [1] 99
## Time difference of 0.2757499 secs
## [1] 100
## Time difference of 0.2756822 secs
## [1] "Build return matrix"
## [1] "Done"
## Time difference of 28.25163 secs## Getting EFS Result in Plotable Format
long = ensemble_result %>%
as.data.frame() %>%
tibble::rownames_to_column("model") %>%
pivot_longer(2:length(.),names_to = "Gene") %>%
#converting genes back to proper notation.
mutate(Gene = gsub("_","-",Gene))
# Ranking Biomarkers
top = long %>%
group_by(Gene) %>%
dplyr::summarise(sum = sum(value)) %>%
dplyr::arrange(-sum)
# Extracting the top 10 ranked biomarkers
topn = top[1:10,] %>%
pull(Gene)
# Filtering the EFS data for plotting purposes
filtered = long %>%
filter(Gene %in% topn)
# # Generating the EFS plot.
# f2d = filtered %>%
# ggplot(aes(reorder(Gene,value),value,fill = model))+
# geom_col()+
# coord_flip()+
# viridis::scale_fill_viridis(discrete = TRUE)+
# ggtitle("Top 10 Performing Protein Biomarkers")+
# xlab("Protein Biomarkers")
#
#
# f2d
#ggsave("02_figures/f2_infection_v_healthy_proteomics/f2d.pdf",f2d,height = 8, width = 8)
# Writting to file to use in supplementary figures.
topn_proteomics = tibble(Gene = topn)
write_csv(topn_proteomics,"01_r_processed_data/top10_infected_v_healthy_proteomics_biomarkers.csv")Plotting
# Reading in data
protein_sum = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
# Processing data
ROC_data = dplyr::filter(protein_sum,Gene %in% topn) %>%
dplyr::mutate(Condition = case_when(Condition %in% c("Faecalis","Faecium") ~ "infected",
TRUE ~ Condition))
source("functions.R")
bottom_plot = make_roc_proteomics(ROC_data,topn[1:2],ncol = 2)
top_plot = make_violin_plot_proteomics(ROC_data,topn[1:2], ncol = 2,pname = "p.signif")
f2c = ggpubr::ggarrange(top_plot,bottom_plot,nrow = 2)
ggsave("02_figures/f2_infection_v_healthy_proteomics/f2c.pdf",f2c,height = 8, width = 8)Figure 3- Metabolomic differences between Healthy and Infected EcB
Processing Metabolomic data
### Unnormalized data ###
#reading and formatting the metabolomics data
data = read_csv("00_r_input_data/mzmine_output/features_quant.csv") %>%
janitor::clean_names() %>%
select(1:3,14:length(.)) %>%
pivot_longer(4:length(.),names_to = "metabolomics_id") %>%
mutate(metabolomics_id = gsub("_mz.*","",metabolomics_id),
metabolomics_id = toupper(metabolomics_id)) %>%
#need to add the RID prefix because having numeric column names (which we will have after pivoting the data) messes up lots of the functions.
mutate(row_id = paste0("RID_",row_id))
# Reading in the mapping
map = read_csv("00_r_input_data/metabolomics_mapping.csv")
md = read_csv("00_r_input_data/clinical_metadata.csv")
# combining the mapping and metadata
map_md = inner_join(map,md,by= "sample_id")
# Combining everything
all = inner_join(data,map_md,by = "metabolomics_id")
write_csv(all,"01_r_processed_data/unnormalized_metabolomics.csv")
##### Normalizing the metabolomics data ######
standard = all %>%
filter(row_id == "RID_6043") %>%
select(metabolomics_id,standard_value = value)
#The normalization that Robert used was to divide the value by the standard value, multiply by 1E6 ( so it isn't negative) and then take the log10 of it.
norm = all %>%
inner_join(standard, by = "metabolomics_id") %>%
mutate(norm_value = log10((value/standard_value) * 1E6))
inf = norm %>%
filter(norm_value == -Inf) %>%
select(row_id,value,standard_value,norm_value) %>%
pull(row_id) %>%
unique()
`%!in%` = Negate(`%in%`)
normalized_data = norm %>%
filter(row_id %!in% inf,
row_id != "RID_6043")
write_csv(normalized_data,"01_r_processed_data/normalized_metabolomics.csv")A - Volcano plots of differentially abundant features.
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv",na = c("N/A","NA"))
## Reading in the normalized metabolomics data ##
normalized_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
# Performing statistics so we can make a volcano plot
stat = normalized_data %>%
dplyr::group_by(row_id) %>%
rstatix::t_test(norm_value ~ infection_status) %>%
rstatix::adjust_pvalue(p.col = "p",output.col = "p.adj_fdr", method = "fdr")
# Extracting the values for each feature for healthy patients
healthy = normalized_data %>%
filter(infection_status == "uninfected") %>%
group_by(row_id) %>%
summarise(healthy = mean(norm_value))
# Extracting the values for each feature for infected patients
infected = normalized_data %>%
filter(infection_status == "infected") %>%
group_by(row_id) %>%
summarise(infected = mean(norm_value))
## Combining the healthy and infected so we can caclulated the log2fc ##
log2fc = inner_join(healthy,infected,by = "row_id") %>%
mutate(log2fc_infected_healthy = log2(infected/healthy))
## Joining the statistics and log2 fold change ##
vp = inner_join(stat,log2fc,by = "row_id") %>%
mutate(annotated = case_when(row_id %in% gnps_annotations$row_id ~ TRUE,
TRUE ~ FALSE))
f3a = vp %>%
ggplot(aes(log2fc_infected_healthy,-log10(p.adj_fdr),color = annotated))+
geom_point(size = 0.5)+
geom_hline(yintercept = -log10(0.05),linetype = "dashed",color = "red")+
xlab("log2FC (Infected/Healthy)")+
ggprism::theme_prism(base_size = 10)+
theme(legend.position = "none")
f3a
ggsave(f3a,filename = "02_figures/f3_infection_v_healthy_metabolomics/f3a.pdf",height = 4, width = 4)B - breakdown on number of features ID
## Reading in the data ##
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv",na = c("N/A","NA"))
normalized_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
unnormalized_data = read_csv("01_r_processed_data/unnormalized_metabolomics.csv")
## Calculating the number of features ##
n_identified_features_gnps = nrow(gnps_annotations)
total_features = length(unique(unnormalized_data$row_id))
features_after_norm = length(unique(normalized_data$row_id))
percent_of_total_features_id = round(n_identified_features_gnps / total_features * 100,)
percent_of_norm_features_id = round(n_identified_features_gnps / features_after_norm * 100,0)
# Putting together a data frame #
feature_count = tibble(`# total features` = total_features,
`# features after normalization` = features_after_norm,
`# of identified features by GNPS` = n_identified_features_gnps,
`% of total features ided` = percent_of_total_features_id,
`% of normalized features ided` = percent_of_norm_features_id) %>%
pivot_longer(1:length(.))
#Making a plot of the table
pdf("02_figures/f3_infection_v_healthy_metabolomics/f3b.pdf",height = 4, width = 4)
gridExtra::grid.table(feature_count)
dev.off()## png
## 2C - enrichment analysis of Metabolites in Infected Samples
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv")
## Structuring the data for enrichment Analysis ##
gnps_enrichment_data = gnps_annotations %>%
#Manual conversion of unreasonable name to what it is listed as in MASSBANK
select(compound_name_cleaned,npclassifier_superclass,npclassifier_class,npclassifier_pathway) %>%
mutate(GO_id = paste0(npclassifier_superclass,"_",npclassifier_class,"_",npclassifier_pathway))
#Making a makeshift "GO" type database using the GNPS annotations.
TERM2GENE = gnps_enrichment_data %>%
dplyr::select(TERM = GO_id,GENE = compound_name_cleaned) %>%
as.data.frame()
TERM2NAME = gnps_enrichment_data %>%
select(-compound_name_cleaned) %>%
pivot_longer(cols = 1:3) %>%
select(keys = GO_id,columns = value) %>%
as.data.frame()
## Extracting the row_ids that were significant from the volcano plot ##
stat_different = vp %>%
filter(p.adj_fdr <= 0.05)
#
#converting RID to compound_name
stat_different_map = inner_join(vp,gnps_annotations,by = "row_id")
## Extracting the row ids that were upregulated in
up = stat_different_map %>%
filter(log2fc_infected_healthy > 0 ) %>%
pull(compound_name_cleaned)
# Performing an enrichment
enriched_in_infected = clusterProfiler::enricher(up,
universe = gnps_enrichment_data$compound_name_cleaned,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME
)
f3c = clusterProfiler::cnetplot(enriched_in_infected,cex_label_category = 1,cex_label_gene = 0.5 )+
theme(legend.position = "none")
f3c
ggsave("02_figures/f3_infection_v_healthy_metabolomics/f3c.pdf",f3c,height = 3,width = 4)down = stat_different_map %>%
filter(log2fc_infected_healthy < 0 ) %>%
pull(compound_name_cleaned)
# Performing an enrichment
enriched_in_healthy = clusterProfiler::enricher(down,
universe = gnps_enrichment_data$compound_name_cleaned,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME,pAdjustMethod = "none"
)
enriched_in_healthy## #
## # over-representation test
## #
## #...@organism UNKNOWN
## #...@ontology UNKNOWN
## #...@gene chr [1:103] "traumatic acid" "methionine" "pentapropylene glycol" ...
## #...pvalues adjusted by 'none' with cutoff <0.05
## #...0 enriched terms found
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141^ There are no enriched terms here.
Staph data for supplement
staph_metabolomics = readxl::read_excel("00_r_input_data/staph_metabolomics.xlsx",skip = 1)
## Structuring the data for enrichment Analysis ##
gnps_enrichment_data = staph_metabolomics %>%
#Manual conversion of unreasonable name to what it is listed as in MASSBANK
select(Unique_ID,Superclass,Class,Subclass) %>%
mutate(GO_id = paste0(Superclass,"_",Class,"_",Subclass))
#Making a makeshift "GO" type database using the GNPS annotations.
TERM2GENE = gnps_enrichment_data %>%
dplyr::select(TERM = GO_id,GENE = Unique_ID) %>%
as.data.frame()
TERM2NAME = gnps_enrichment_data %>%
select(-Unique_ID) %>%
pivot_longer(cols = 1:3) %>%
select(keys = GO_id,columns = value) %>%
as.data.frame()
## Extracting the row_ids that were significant from the volcano plot ##
stat_different = staph_metabolomics %>%
mutate(Infection_pvalue = as.numeric(Infection_pvalue)) %>%
filter(Infection_pvalue <= 0.05)
## Extracting the row ids that were upregulated in
up = stat_different %>%
filter(Infection_FC > 1) %>%
pull(Unique_ID)
# Performing an enrichment
enriched_in_infected = clusterProfiler::enricher(up,
universe = gnps_enrichment_data$Unique_ID,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME
)
## Here is the data frame with only the identified unique IDs
names = staph_metabolomics %>%
filter(Unique_ID %in% c(2510,2915,1917,2800,3125,3138,2263,2240,2470,3507,2767,1898,3509,1944))
# Manually converting the RID to reasonable names
fixed_names_vector = c("2510",
"12-Ketodeoxycholic acid",
"2-((4R)-4-((3R,5S,7R,9S,10S,13R,14S,17R)-3,7-dihydroxy-10,13-dimethylhexadecahydro-1H-cyclopenta[a]phenanthren-17-yl)pentanamido)ethane-1-sulfonic acid",
"2-((4R)-4-((3R,5S,7R,9S,10S,13R,14S,17R)-3,7-dihydroxy-10,13-dimethylhexadecahydro-1H-cyclopenta[a]phenanthren-17-yl)pentanamido)ethane-1-sulfonic acid",
"(2S,6R)-2-methyl-6-((3R,5S,7R,9S,10S,12S,13R,14S,17R)-3,7,12-trihydroxy-10,13-dimethylhexadecahydro-1H-cyclopenta[a]phenanthren-17-yl)heptanoic acid",
"3138",
"2-((4R)-4-((3R,5S,7R,9S,10S,13R,14S,17R)-3,7-dihydroxy-10,13-dimethylhexadecahydro-1H-cyclopenta[a]phenanthren-17-yl)pentanamido)ethane-1-sulfonic acid",
"Tauroursodeoxycholic acid",
"3beta-Hydroxy-5-cholenoic acid",
"Taurocholic acid",
"2767",
"Taurocholic acid",
"Taurocholic acid",
"Glycocholic acid")
fixed_names = paste0(fixed_names_vector,collapse = "/")
enriched_in_infected@result$geneID = fixed_names
up_plot = clusterProfiler::cnetplot(enriched_in_infected,cex_label_category = 1, cex_label_gene = 0.4)+
theme(legend.position = "none")
up_plot
ggsave("03_supplemental/s2_networks_compared_to_healthy/up_in_staph_metabolomics.pdf",up_plot,width = 4, height = 4)## Extracting the row ids that were upregulated in
down = stat_different %>%
filter(Infection_FC < 1) %>%
pull(Unique_ID)
# Performing an enrichment\
# Performing an enrichment
enriched_in_healthy = clusterProfiler::enricher(down,
universe = gnps_enrichment_data$Unique_ID,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME
)
## Here is the data frame with only the identified unique IDs
names = staph_metabolomics %>%
filter(Unique_ID %in% c(554,132,450,1939,33,372,99,557,1057,586,434,128))
fixed_names_vector = c("1-Myristoyl-sn-glycero-3-phosphocholine",
"132",
"1-Octadecyl-2-acetyl-sn-glycero-3-phosphocholine",
"1-Myristoyl-sn-glycero-3-phosphocholine",
"1-Pentadecanoyl-sn-glycero-3-phosphocholine",
"372",
"1-Octadecyl-sn-glycero-3-phosphocholine",
"1,2-Dioctanoyl PC",
"1-(1Z-Octadecenyl)-sn-glycero-3-phosphocholine",
"1-(1Z-Octadecenyl)-sn-glycero-3-phosphocholine",
"434",
"1-O-Hexadecyl-2-O-(2E-butenoyl)-sn-glyceryl-3-phosphocholine")
fixed_names = paste0(fixed_names_vector,collapse = "/")
enriched_in_healthy@result$geneID = fixed_names
healthy_plot = clusterProfiler::cnetplot(enriched_in_healthy,cex_label_category = 1,show_category = 100,cex_label_gene = 0.5)+
theme(legend.position = "none")
healthy_plot
ggsave("03_supplemental/s2_networks_compared_to_healthy/down_in_staph_metabolomics.pdf",healthy_plot,width = 4, height = 4)D - Assesment of top biomarkers
Data processing
## Reading in the normalized data ##
normalized_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv")
## Switching to the wide format ###
norm_wide = normalized_data %>%
select(-row_m_z,-row_retention_time,-value) %>%
tidyr::pivot_wider(names_from = row_id,values_from = norm_value)
## Looking for only columns that were annotated ##
gnps_and_in_norm = intersect(unique(gnps_annotations$row_id), colnames(norm_wide))
gnps_annotated_metabolites = norm_wide %>%
as.data.frame() %>%
select(condition,gnps_and_in_norm)
## Extracting the numeric data ##
efs_sel = gnps_annotated_metabolites %>%
mutate(classlabels = case_when(condition %in% c("faecalis","faecium") ~ 0,
TRUE ~ 1)) %>%
select(classlabels,2:length(.)) %>%
#absolutely needs to be a data frame and not a tibble for EFS package
as.data.frame()
# Seting seed incase this isn't done internally in EFS
set.seed(2)
## Running EFS feature selection ##
ensemble_result = EFS::ensemble_fs(efs_sel,
classnumber = 1,cor_threshold = 0)## [1] "default value for NA_threshold = 0.2"
## [1] "default value for runs = 100"
## [1] "default value for selection is c(TRUE, TRUE, TRUE,TRUE, TRUE, TRUE, FALSE, FALSE)"
## [1] "Start Median"
## [1] "Start Pearson"
## [1] "Start Spearman"
## [1] "Start LogReg"
## [1] "Start RF"
## [1] 1
## Time difference of 0.2077355 secs
## [1] 2
## Time difference of 0.2080009 secs
## [1] 3
## Time difference of 0.2202313 secs
## [1] 4
## Time difference of 0.2076793 secs
## [1] 5
## Time difference of 0.2077954 secs
## [1] 6
## Time difference of 0.2077839 secs
## [1] 7
## Time difference of 0.2060838 secs
## [1] 8
## Time difference of 0.2060761 secs
## [1] 9
## Time difference of 0.2079167 secs
## [1] 10
## Time difference of 0.2071259 secs
## [1] 11
## Time difference of 0.2068918 secs
## [1] 12
## Time difference of 0.206836 secs
## [1] 13
## Time difference of 0.2063446 secs
## [1] 14
## Time difference of 0.205981 secs
## [1] 15
## Time difference of 0.2066832 secs
## [1] 16
## Time difference of 0.2066953 secs
## [1] 17
## Time difference of 0.204138 secs
## [1] 18
## Time difference of 0.2049105 secs
## [1] 19
## Time difference of 0.2200446 secs
## [1] 20
## Time difference of 0.2067101 secs
## [1] 21
## Time difference of 0.2070434 secs
## [1] 22
## Time difference of 0.2078326 secs
## [1] 23
## Time difference of 0.2070866 secs
## [1] 24
## Time difference of 0.2069833 secs
## [1] 25
## Time difference of 0.2058172 secs
## [1] 26
## Time difference of 0.2066183 secs
## [1] 27
## Time difference of 0.206538 secs
## [1] 28
## Time difference of 0.2069356 secs
## [1] 29
## Time difference of 0.2066898 secs
## [1] 30
## Time difference of 0.2053015 secs
## [1] 31
## Time difference of 0.2070882 secs
## [1] 32
## Time difference of 0.2069449 secs
## [1] 33
## Time difference of 0.2059269 secs
## [1] 34
## Time difference of 0.219466 secs
## [1] 35
## Time difference of 0.2070279 secs
## [1] 36
## Time difference of 0.2060678 secs
## [1] 37
## Time difference of 0.2063031 secs
## [1] 38
## Time difference of 0.2069743 secs
## [1] 39
## Time difference of 0.2070875 secs
## [1] 40
## Time difference of 0.206183 secs
## [1] 41
## Time difference of 0.2072814 secs
## [1] 42
## Time difference of 0.2084429 secs
## [1] 43
## Time difference of 0.2059329 secs
## [1] 44
## Time difference of 0.2053804 secs
## [1] 45
## Time difference of 0.2084107 secs
## [1] 46
## Time difference of 0.2081971 secs
## [1] 47
## Time difference of 0.2062514 secs
## [1] 48
## Time difference of 0.2064102 secs
## [1] 49
## Time difference of 0.2201416 secs
## [1] 50
## Time difference of 0.2057588 secs
## [1] 51
## Time difference of 0.2058947 secs
## [1] 52
## Time difference of 0.2073598 secs
## [1] 53
## Time difference of 0.2071331 secs
## [1] 54
## Time difference of 0.2076197 secs
## [1] 55
## Time difference of 0.2073579 secs
## [1] 56
## Time difference of 0.2066422 secs
## [1] 57
## Time difference of 0.2079682 secs
## [1] 58
## Time difference of 0.2065983 secs
## [1] 59
## Time difference of 0.2045555 secs
## [1] 60
## Time difference of 0.2074871 secs
## [1] 61
## Time difference of 0.2073925 secs
## [1] 62
## Time difference of 0.2081804 secs
## [1] 63
## Time difference of 0.2078836 secs
## [1] 64
## Time difference of 0.2189724 secs
## [1] 65
## Time difference of 0.2077644 secs
## [1] 66
## Time difference of 0.2067437 secs
## [1] 67
## Time difference of 0.2084954 secs
## [1] 68
## Time difference of 0.2074072 secs
## [1] 69
## Time difference of 0.2073555 secs
## [1] 70
## Time difference of 0.2065432 secs
## [1] 71
## Time difference of 0.2057145 secs
## [1] 72
## Time difference of 0.2071922 secs
## [1] 73
## Time difference of 0.2075787 secs
## [1] 74
## Time difference of 0.2082648 secs
## [1] 75
## Time difference of 0.2077382 secs
## [1] 76
## Time difference of 0.2053576 secs
## [1] 77
## Time difference of 0.2068615 secs
## [1] 78
## Time difference of 0.2075274 secs
## [1] 79
## Time difference of 0.2205014 secs
## [1] 80
## Time difference of 0.2076635 secs
## [1] 81
## Time difference of 0.2067199 secs
## [1] 82
## Time difference of 0.2069244 secs
## [1] 83
## Time difference of 0.2076631 secs
## [1] 84
## Time difference of 0.2075896 secs
## [1] 85
## Time difference of 0.2066407 secs
## [1] 86
## Time difference of 0.2066698 secs
## [1] 87
## Time difference of 0.2065969 secs
## [1] 88
## Time difference of 0.2061479 secs
## [1] 89
## Time difference of 0.205853 secs
## [1] 90
## Time difference of 0.2073526 secs
## [1] 91
## Time difference of 0.2053549 secs
## [1] 92
## Time difference of 0.2048538 secs
## [1] 93
## Time difference of 0.2067351 secs
## [1] 94
## Time difference of 0.2182655 secs
## [1] 95
## Time difference of 0.2076652 secs
## [1] 96
## Time difference of 0.2089598 secs
## [1] 97
## Time difference of 0.2074554 secs
## [1] 98
## Time difference of 0.2068152 secs
## [1] 99
## Time difference of 0.2063401 secs
## [1] 100
## Time difference of 0.2064116 secs
## [1] "Build return matrix"
## [1] "Done"
## Time difference of 20.90489 secs## Formating the data as long ##
long = ensemble_result %>%
as.data.frame() %>%
tibble::rownames_to_column("model") %>%
pivot_longer(2:length(.),names_to = "row_id")
## Extracting top metabolites ##
top = long %>%
group_by(row_id) %>%
dplyr::summarise(sum = sum(value)) %>%
dplyr::arrange(-sum)
## Pulling top N metabolites ##
topn = top[1:10,] %>%
pull(row_id)
#writing file contining the topn row ids
topn_metabolomics = tibble(row_id = topn)
write_csv(topn_metabolomics,"01_r_processed_data/top10_infected_v_healthy_metabolomics_biomarkers.csv")
## Filtering the data, combining annotation names ##
filtered = long %>%
filter(row_id %in% topn) %>%
inner_join(gnps_annotations,by = "row_id")
## Plotting the data ##
# f3d = filtered %>%
# ggplot(aes(reorder(compound_name_cleaned,value),value,fill = model))+
# geom_col()+
# coord_flip()+
# viridis::scale_fill_viridis(discrete = TRUE)+
# ggtitle("Top 10 Performing Metabolite Biomarkers of Infection")+
# xlab("Metabolite Biomarkers")
#
# f3d
#ggsave("02_figures/f3_infection_v_healthy_metabolomics/f3d.pdf",f3d,height = 8,width = 8,units = "in")Plotting
## Arranging Data so it is in the correct order ###
ROC_data = dplyr::filter(normalized_data,row_id %in% topn[1:2]) %>%
dplyr::arrange(factor(row_id,topn)) %>%
inner_join(gnps_annotations,by = c("row_id" = "row_id")) %>%
mutate(condition = case_when(condition == "healthy" ~ "healthy",
TRUE ~ "infected"))
source("functions.R")
top = make_violin_plot_metabolomics(ROC_data,topn = topn[1:2],ncol = 2,pname = "p.signif")
bottom = make_roc_metabolomics(ROC_data,ncol = 2)
f3d = ggpubr::ggarrange(top,bottom,nrow = 2)
f3d
ggsave("02_figures/f3_infection_v_healthy_metabolomics/f3d.pdf",f3d,height=8,width=8)Figure 4- Proteomic Differences between Enterococcus Faecium and Faecalis
A- Volcano plot
faecalis_v_faecium = readr::read_csv("01_r_processed_data/faecalis_v_faecium.csv") %>%
left_join(protein_to_gene,by = c("uniprotid" = "ProteinID"))
f4a = faecalis_v_faecium %>%
ggplot(aes(log2FC,-log10(adj.pvalue)))+
geom_point(size = 0.5)+
geom_hline(yintercept = -log10(0.05),linetype = "dashed",color = "red")+
xlab("Log2FC Faecalis / Faecium")+
ggprism::theme_prism(base_size = 10)
f4a
ggsave("02_figures/f4_faecalis_v_faecium_proteomics/f4a.pdf",f4a,height = 4, width = 4)B Network analysis up in Faecalis
up_in_faecalis = faecalis_v_faecium %>%
filter(log2FC>0) %>%
filter(adj.pvalue <= 0.05) %>%
arrange(adj.pvalue) %>%
pull(Gene)
#Getting all the go annotations
go = read_csv("00_r_input_data/go_annotations.csv") %>%
dplyr::select(1,length(.)) %>%
tidyr::separate_rows(Gene.Ontology.IDs, sep = ';') %>%
dplyr::mutate(Gene.Ontology.IDs = gsub(" ","",Gene.Ontology.IDs)) %>%
# Going from Protein ID to gene ID
left_join(protein_to_gene,by = c("Entry" = "ProteinID")) %>%
select(-Entry) %>%
select(Entry = Gene,Gene.Ontology.IDs)
#usng enricher
TERM2GENE = go %>%
dplyr::select(TERM = Gene.Ontology.IDs,GENE = Entry) %>%
as.data.frame()
TERM2NAME = AnnotationDbi::select(GO.db::GO.db,keys = unique(TERM2GENE$TERM),columns = c("TERM")) %>%
as.data.frame()
enrich = clusterProfiler::enricher(up_in_faecalis,
universe = go$Entry,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME)
f4b = clusterProfiler::cnetplot(enrich,
categorySize="pvalue",
cex_label_category = 1,
cex_label_gene = 0.5,
showCategory = 100)+
theme(legend.position = "none")
f4b
ggsave("02_figures/f4_faecalis_v_faecium_proteomics/f4b.pdf",height = 4, width =8)C - Immunoglobulin comparisons
protein_sum = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
anot = read_csv("00_r_input_data/all_protein_annotations.csv")
sum_annotated = left_join(protein_sum,anot,by = c("uniprotid" = "Entry"))
ig = sum_annotated %>%
filter(grepl(".*Immunoglobulin.*",Protein.names)) %>%
filter(!grepl(".*receptor.*",Protein.names)) %>%
#Removing one problematic immunglobulin. Negative values for 3 of the samples messed things up.
#Doesn't change the story, just the plot.
filter(uniprotid != "A0A0B4J2B5")
stat = ig %>%
rstatix::t_test(Abundance ~ Condition)
f4c = ig %>%
ggplot(aes(Condition,Abundance),size = 0.5)+
geom_violin()+
geom_boxplot()+
ggpubr::stat_pvalue_manual(stat,
y.position = 22,
step.increase = 0.1)+
ggprism::theme_prism(base_size = 10)
f4c
ggsave("02_figures/f4_faecalis_v_faecium_proteomics/f4c.pdf",f4c,height = 3.8,width = 4)D - terms up in faecium.
up_in_faecium = faecalis_v_faecium %>%
filter(log2FC<0) %>%
filter(adj.pvalue <= 0.05) %>%
arrange(adj.pvalue) %>%
pull(Gene)
enrich = clusterProfiler::enricher(up_in_faecium,
universe = go$Entry,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME)
f4d = clusterProfiler::cnetplot(enrich,
categorySize="pvalue",
cex_label_category = 1,
cex_label_gene = 0.5,
showCategory = 100)+
theme(legend.position = "none")
f4d
ggsave("02_figures/f4_faecalis_v_faecium_proteomics/f4d.pdf",f4d,height = 3.8, width = 4)E - Assesment of Biomarkers
Data Processing
library(EFS)
pl = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
clean = pl %>%
dplyr::group_by(Mixture,Gene,Channel,Condition) %>%
dplyr::summarise(sum = sum(Abundance))
faecium_v_faecalis_l = clean %>%
filter(Condition %in% c("Faecium", "Faecalis")) %>%
mutate(classlabels = case_when(Condition == "Faecalis" ~ 1,
TRUE ~ 0)) %>%
# EFS can't handle the - charachter.
mutate(Gene = gsub("-","_",Gene))
faecium_v_faecalis_w = faecium_v_faecalis_l %>%
pivot_wider(names_from = Gene,
values_from = sum) %>%
mutate(sampleid = paste0(Mixture,".",Channel)) %>%
ungroup() %>%
dplyr::select(-Mixture,-Channel,-Condition,-sampleid) %>%
select_if(~ !any(is.na(.))) %>%
as.data.frame()
#need to come up with some sort of way to combine all of the fractions and plexes
# The problem is that na.omit filters everything out.
ensemble_result = EFS::ensemble_fs(faecium_v_faecalis_w,1,cor_threshold = 0)## [1] "default value for NA_threshold = 0.2"
## [1] "default value for runs = 100"
## [1] "default value for selection is c(TRUE, TRUE, TRUE,TRUE, TRUE, TRUE, FALSE, FALSE)"
## [1] "Start Median"
## [1] "Start Pearson"
## [1] "Start Spearman"
## [1] "Start LogReg"
## [1] "Start RF"
## [1] 1
## Time difference of 0.3044047 secs
## [1] 2
## Time difference of 0.3036761 secs
## [1] 3
## Time difference of 0.3017886 secs
## [1] 4
## Time difference of 0.3019941 secs
## [1] 5
## Time difference of 0.3029802 secs
## [1] 6
## Time difference of 0.3162594 secs
## [1] 7
## Time difference of 0.302953 secs
## [1] 8
## Time difference of 0.3029811 secs
## [1] 9
## Time difference of 0.3030562 secs
## [1] 10
## Time difference of 0.3023045 secs
## [1] 11
## Time difference of 0.3040721 secs
## [1] 12
## Time difference of 0.3031714 secs
## [1] 13
## Time difference of 0.3033929 secs
## [1] 14
## Time difference of 0.3036571 secs
## [1] 15
## Time difference of 0.3033013 secs
## [1] 16
## Time difference of 0.3042903 secs
## [1] 17
## Time difference of 0.302063 secs
## [1] 18
## Time difference of 0.3011875 secs
## [1] 19
## Time difference of 0.3035676 secs
## [1] 20
## Time difference of 0.3045449 secs
## [1] 21
## Time difference of 0.3038983 secs
## [1] 22
## Time difference of 0.3061092 secs
## [1] 23
## Time difference of 0.305038 secs
## [1] 24
## Time difference of 0.306289 secs
## [1] 25
## Time difference of 0.306392 secs
## [1] 26
## Time difference of 0.3044379 secs
## [1] 27
## Time difference of 0.3163276 secs
## [1] 28
## Time difference of 0.3042281 secs
## [1] 29
## Time difference of 0.3053284 secs
## [1] 30
## Time difference of 0.3038881 secs
## [1] 31
## Time difference of 0.3043685 secs
## [1] 32
## Time difference of 0.3034108 secs
## [1] 33
## Time difference of 0.3032651 secs
## [1] 34
## Time difference of 0.3030298 secs
## [1] 35
## Time difference of 0.3044555 secs
## [1] 36
## Time difference of 0.3043308 secs
## [1] 37
## Time difference of 0.3042214 secs
## [1] 38
## Time difference of 0.3045068 secs
## [1] 39
## Time difference of 0.3024266 secs
## [1] 40
## Time difference of 0.3037584 secs
## [1] 41
## Time difference of 0.3049035 secs
## [1] 42
## Time difference of 0.3033485 secs
## [1] 43
## Time difference of 0.3019032 secs
## [1] 44
## Time difference of 0.3026931 secs
## [1] 45
## Time difference of 0.3034542 secs
## [1] 46
## Time difference of 0.3043036 secs
## [1] 47
## Time difference of 0.3040738 secs
## [1] 48
## Time difference of 0.3168948 secs
## [1] 49
## Time difference of 0.3042774 secs
## [1] 50
## Time difference of 0.3038578 secs
## [1] 51
## Time difference of 0.3030453 secs
## [1] 52
## Time difference of 0.3036082 secs
## [1] 53
## Time difference of 0.3027771 secs
## [1] 54
## Time difference of 0.3045697 secs
## [1] 55
## Time difference of 0.302824 secs
## [1] 56
## Time difference of 0.3025713 secs
## [1] 57
## Time difference of 0.3026772 secs
## [1] 58
## Time difference of 0.3041263 secs
## [1] 59
## Time difference of 0.302392 secs
## [1] 60
## Time difference of 0.3045285 secs
## [1] 61
## Time difference of 0.3017957 secs
## [1] 62
## Time difference of 0.3033764 secs
## [1] 63
## Time difference of 0.3044672 secs
## [1] 64
## Time difference of 0.3029482 secs
## [1] 65
## Time difference of 0.3042769 secs
## [1] 66
## Time difference of 0.3030345 secs
## [1] 67
## Time difference of 0.3026719 secs
## [1] 68
## Time difference of 0.3057072 secs
## [1] 69
## Time difference of 0.3158166 secs
## [1] 70
## Time difference of 0.3043897 secs
## [1] 71
## Time difference of 0.3044121 secs
## [1] 72
## Time difference of 0.3027127 secs
## [1] 73
## Time difference of 0.3043251 secs
## [1] 74
## Time difference of 0.3041878 secs
## [1] 75
## Time difference of 0.303539 secs
## [1] 76
## Time difference of 0.3031757 secs
## [1] 77
## Time difference of 0.301722 secs
## [1] 78
## Time difference of 0.3041215 secs
## [1] 79
## Time difference of 0.3041501 secs
## [1] 80
## Time difference of 0.3034911 secs
## [1] 81
## Time difference of 0.3033597 secs
## [1] 82
## Time difference of 0.3028359 secs
## [1] 83
## Time difference of 0.3023705 secs
## [1] 84
## Time difference of 0.3025551 secs
## [1] 85
## Time difference of 0.3020325 secs
## [1] 86
## Time difference of 0.3040209 secs
## [1] 87
## Time difference of 0.3040607 secs
## [1] 88
## Time difference of 0.3040473 secs
## [1] 89
## Time difference of 0.3026786 secs
## [1] 90
## Time difference of 0.3314621 secs
## [1] 91
## Time difference of 0.3026755 secs
## [1] 92
## Time difference of 0.3033009 secs
## [1] 93
## Time difference of 0.3056602 secs
## [1] 94
## Time difference of 0.3056369 secs
## [1] 95
## Time difference of 0.30424 secs
## [1] 96
## Time difference of 0.3032782 secs
## [1] 97
## Time difference of 0.3038797 secs
## [1] 98
## Time difference of 0.3042393 secs
## [1] 99
## Time difference of 0.3024821 secs
## [1] 100
## Time difference of 0.3033736 secs
## [1] "Build return matrix"
## [1] "Done"
## Time difference of 30.55866 secslong = ensemble_result %>%
as.data.frame() %>%
#Converting gene names back
tibble::rownames_to_column("model") %>%
pivot_longer(2:length(.),names_to = "Gene") %>%
mutate(Gene = gsub("_","-",Gene))
top = long %>%
group_by(Gene) %>%
dplyr::summarise(sum = sum(value)) %>%
dplyr::arrange(-sum)
topn = top[1:10,] %>%
pull(Gene)
# writing to file
topn_proteomics = tibble(Gene = topn)
write_csv(topn_proteomics,"01_r_processed_data/top10_faecalis_v_faecium_proteomics_biomarkers.csv")
# f4d = long %>%
# filter(Gene %in% topn) %>%
# ggplot(aes(reorder(Gene,value),value,fill = model))+
# geom_col()+
# coord_flip()+
# viridis::scale_fill_viridis(discrete = TRUE)+
# ggtitle("Top 10 Biomarkers Differentiating Faecalis from Faecium by EFS")+
# xlab("Protein Biomarkers")
#
# f4d
# ggsave("02_figures/f4_faecalis_v_faecium_proteomics/f4d.pdf",f4d,height = 8,width = 8)Plotting
ROC_data = dplyr::filter(protein_sum,Gene %in% topn[1:2]) %>%
dplyr::filter(Condition %in% c("Faecalis","Faecium"))
source("functions.R")
f4d = make_violin_plot_proteomics(ROC_data,topn[1:2],ncol = 2,pname = "p.signif")
f4e = make_roc_proteomics(ROC_data,topn[1:2],ncol=2)
f4d = ggpubr::ggarrange(f4d,f4e,nrow = 2)
ggsave("02_figures/f4_faecalis_v_faecium_proteomics/f4d.pdf",f4d,height=8,width=8)? Are proteins involved in the gene ontology results confounded?
# Extracting the GO terms in
enriched_in_faecium = enrich@result %>%
filter(p.adjust <= 0.05)
#Looking at the proteins in the cholesterol portions of the network
proteins_in_network = c("APOC3","APOE","APOC1","APOC2","PCYOX1","CLU")
clinical_md = read_csv("00_r_input_data/clinical_metadata.csv")
protein_sum = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
# Defining the input metadata. Need to chose which columns to keep here.
md_assesment_data= protein_sum %>%
filter(Condition != "Healthy") %>%
inner_join(clinical_md, by = c("BioReplicate" = "sample_id")) %>%
select(sample_id,Gene,bmi,bacteremia_duration,charleson_index,congestive_heart_failure,
death_during_admission,history_of_myocardial_infarction,transplant_type,pathogen,icu_admission,
hypotension,gender,race,smoking,polymicrobial_bacteremia,neutropenia_anc_1500,leukemia_or_lymphoma,Abundance) %>%
ungroup()
source("functions.R")
stat = metadata_assesment(md_assesment_data,
metric = Gene,
sample_id = sample_id,
features = proteins_in_network,
outcome = Abundance,
outcome_c = "Abundance")
stat_final = rstatix::adjust_pvalue(stat,"p",output.col = "p.adj",method = "fdr")
## Plotting ###
t = stat_final %>%
ggplot(aes(reorder(variable,p.adj),-log10(p.adj)))+
geom_col()+
geom_hline(yintercept = -log10(0.05),linetype = "dashed", color = "red")+
coord_flip()+
facet_wrap(~Gene,nrow = 2) +
xlab("metadata_field")
t
plotlist = list()
for(i in proteins_in_network){
f = filter(md_assesment_data,Gene %in% i,
pathogen %in% c("faecalis","faecium"),
transplant_type == "none")
stat = f %>%
rstatix::t_test(Abundance ~ pathogen,) %>%
#rstatix::add_significance() %>%
rstatix::add_xy_position(fun = "median")
p1 = ggplot(f,aes(pathogen,Abundance)) +
geom_violin()+
geom_point()+
ggtitle(i)+
ggpubr::stat_pvalue_manual(stat,)
plotlist[[i]] = p1
}
library(gridExtra)
t2 = do.call("grid.arrange", plotlist,)
##? Are differentially abundant proteins the same ones that are different when comparing healthy to infected?
healthy_v_infected = readr::read_csv("01_r_processed_data/Infected_vs_Healthy.csv")
up_in_infected = healthy_v_infected %>%
filter(adj.pvalue <= 0.05,
log2FC>0) %>%
pull(uniprotid)
down_in_infected = healthy_v_infected %>%
filter(adj.pvalue <= 0.05,
log2FC<0) %>%
pull(uniprotid)
faecalis_v_faecium = readr::read_csv("01_r_processed_data/faecalis_v_faecium.csv") %>%
left_join(protein_to_gene,by = c("uniprotid" = "ProteinID")) %>%
mutate(also_in = case_when(uniprotid %in% up_in_infected ~ "infected",
uniprotid %in% down_in_infected ~ "healthy",
TRUE ~ "distinct"))
vp = faecalis_v_faecium %>%
ggplot(aes(log2FC,-log10(adj.pvalue),color = also_in))+
geom_point(size = 0.5)+
geom_hline(yintercept = -log10(0.05),linetype = "dashed",color = "red")+
xlab("Log2FC Faecalis / Faecium")+
ggprism::theme_prism(base_size = 10)
vp
Figure 5 - Metabolomic Differences between Faecalis and Faecium Infected samples
A - Volcano plot of data.
## Reading in Data
normalized_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv")
# Performing statistics so we can make a volcano plot
stat = normalized_data %>%
filter(condition %in% c("faecalis","faecium")) %>%
dplyr::group_by(row_id) %>%
rstatix::t_test(norm_value ~ condition) %>%
rstatix::adjust_pvalue(p.col = "p",output.col = "p.adj_fdr", method = "fdr")
# Calculating log2fc
log2fc_data = normalized_data %>%
dplyr::select(metabolomics_id,row_id,condition,norm_value)
# Extracting just the faecalis data
faecalis = log2fc_data %>%
filter(condition == "faecalis") %>%
group_by(row_id) %>%
summarise(faecalis = mean(norm_value))
# Extracting just the faecium data
faecium = log2fc_data %>%
filter(condition == "faecium") %>%
group_by(row_id) %>%
summarise(faecium = mean(norm_value))
# Calculating log2fc
log2fc = inner_join(faecalis,faecium,by = "row_id") %>%
mutate(log2fc_faecalis_faecium = log2(faecalis/faecium))
# Combining log2fc and statistics
vp = inner_join(stat,log2fc,by = "row_id") %>%
mutate(annotated = case_when(row_id %in% gnps_annotations$row_id ~ TRUE,
TRUE ~ FALSE))
annotations = vp %>%
filter(p.adj_fdr <= 0.05,
annotated == TRUE) %>%
inner_join(gnps_annotations, by = "row_id")
f5a = vp %>%
ggplot(aes(log2fc_faecalis_faecium,-log10(p.adj_fdr),color = annotated))+
geom_point(size = 0.5)+
geom_hline(yintercept = -log10(0.05),linetype = "dashed",color = "red")+
xlab("log2FC (Faecalis/Faecium)")+
ggrepel::geom_label_repel(data = annotations,aes(label = compound_name_cleaned,x= log2fc_faecalis_faecium, y = -log10(p.adj_fdr)),size = 2.5,label.size = 0.1)+
ggprism::theme_prism(base_size = 10)+
theme(legend.position = "none")
f5a
ggsave("02_figures/f5_faecalis_v_faecium_metabolomics/f5a.pdf",f5a,height = 3.8, width =4)B - Biomarker Assesment
Data Processing
## Reading in the normalized data ##
normalized_data = read_csv("01_r_processed_data/normalized_metabolomics.csv") %>%
#Removing reduntant RIDs. These were found to be the same as other top biomarkers
filter(row_id != "RID_8856",
row_id != "RID_22529",
row_id != "RID_22270",
row_id != "RID_9436")
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv")
## Switching to the wide format ###
norm_wide = normalized_data %>%
select(-row_m_z,-row_retention_time,-value) %>%
tidyr::pivot_wider(names_from = row_id,values_from = norm_value)
## Looking for only columns that were annotated ##
gnps_and_in_norm = intersect(unique(gnps_annotations$row_id), colnames(norm_wide))
gnps_annotated_metabolites = norm_wide %>%
as.data.frame() %>%
select(condition,gnps_and_in_norm)
## Extracting the numeric data ##
efs_sel = gnps_annotated_metabolites %>%
filter(condition %in% c("faecalis","faecium")) %>%
mutate(classlabels = case_when(condition == "faecalis" ~ 0,
TRUE ~ 1)) %>%
select(classlabels,2:length(.)) %>%
#absolutely needs to be a data frame and not a tibble for EFS package
as.data.frame()
# Seting seed incase this isn't done internally in EFS
set.seed(2)
## Running EFS feature selection ##
ensemble_result = EFS::ensemble_fs(efs_sel,
classnumber = 1,
cor_threshold = 0)## [1] "default value for NA_threshold = 0.2"
## [1] "default value for runs = 100"
## [1] "default value for selection is c(TRUE, TRUE, TRUE,TRUE, TRUE, TRUE, FALSE, FALSE)"
## [1] "Start Median"
## [1] "Start Pearson"
## [1] "Start Spearman"
## [1] "Start LogReg"
## [1] "Start RF"
## [1] 1
## Time difference of 0.2516742 secs
## [1] 2
## Time difference of 0.2498114 secs
## [1] 3
## Time difference of 0.2495921 secs
## [1] 4
## Time difference of 0.2485569 secs
## [1] 5
## Time difference of 0.2488589 secs
## [1] 6
## Time difference of 0.2495742 secs
## [1] 7
## Time difference of 0.2488408 secs
## [1] 8
## Time difference of 0.2488775 secs
## [1] 9
## Time difference of 0.2494268 secs
## [1] 10
## Time difference of 0.2502444 secs
## [1] 11
## Time difference of 0.2486207 secs
## [1] 12
## Time difference of 0.2500749 secs
## [1] 13
## Time difference of 0.2486598 secs
## [1] 14
## Time difference of 0.2491944 secs
## [1] 15
## Time difference of 0.2492626 secs
## [1] 16
## Time difference of 0.2567317 secs
## [1] 17
## Time difference of 0.2493773 secs
## [1] 18
## Time difference of 0.2493172 secs
## [1] 19
## Time difference of 0.2486112 secs
## [1] 20
## Time difference of 0.2494423 secs
## [1] 21
## Time difference of 0.2490165 secs
## [1] 22
## Time difference of 0.2487645 secs
## [1] 23
## Time difference of 0.250025 secs
## [1] 24
## Time difference of 0.2488635 secs
## [1] 25
## Time difference of 0.2486036 secs
## [1] 26
## Time difference of 0.2489514 secs
## [1] 27
## Time difference of 0.2481308 secs
## [1] 28
## Time difference of 0.2489939 secs
## [1] 29
## Time difference of 0.2487392 secs
## [1] 30
## Time difference of 0.2478318 secs
## [1] 31
## Time difference of 0.2546539 secs
## [1] 32
## Time difference of 0.2493443 secs
## [1] 33
## Time difference of 0.2488678 secs
## [1] 34
## Time difference of 0.2497747 secs
## [1] 35
## Time difference of 0.2498782 secs
## [1] 36
## Time difference of 0.2485714 secs
## [1] 37
## Time difference of 0.2490556 secs
## [1] 38
## Time difference of 0.2486086 secs
## [1] 39
## Time difference of 0.2490003 secs
## [1] 40
## Time difference of 0.2470076 secs
## [1] 41
## Time difference of 0.2482014 secs
## [1] 42
## Time difference of 0.2477937 secs
## [1] 43
## Time difference of 0.2478981 secs
## [1] 44
## Time difference of 0.2480958 secs
## [1] 45
## Time difference of 0.248702 secs
## [1] 46
## Time difference of 0.24968 secs
## [1] 47
## Time difference of 0.2551548 secs
## [1] 48
## Time difference of 0.2492054 secs
## [1] 49
## Time difference of 0.2492383 secs
## [1] 50
## Time difference of 0.2499428 secs
## [1] 51
## Time difference of 0.2482898 secs
## [1] 52
## Time difference of 0.2483733 secs
## [1] 53
## Time difference of 0.2490816 secs
## [1] 54
## Time difference of 0.2495644 secs
## [1] 55
## Time difference of 0.2493949 secs
## [1] 56
## Time difference of 0.2482359 secs
## [1] 57
## Time difference of 0.2495549 secs
## [1] 58
## Time difference of 0.2490742 secs
## [1] 59
## Time difference of 0.2481968 secs
## [1] 60
## Time difference of 0.2487226 secs
## [1] 61
## Time difference of 0.2486165 secs
## [1] 62
## Time difference of 0.2491362 secs
## [1] 63
## Time difference of 0.2539854 secs
## [1] 64
## Time difference of 0.2487485 secs
## [1] 65
## Time difference of 0.2487135 secs
## [1] 66
## Time difference of 0.249769 secs
## [1] 67
## Time difference of 0.2487483 secs
## [1] 68
## Time difference of 0.2499466 secs
## [1] 69
## Time difference of 0.2482245 secs
## [1] 70
## Time difference of 0.2493658 secs
## [1] 71
## Time difference of 0.2497916 secs
## [1] 72
## Time difference of 0.248513 secs
## [1] 73
## Time difference of 0.2484546 secs
## [1] 74
## Time difference of 0.2488077 secs
## [1] 75
## Time difference of 0.2490096 secs
## [1] 76
## Time difference of 0.2485814 secs
## [1] 77
## Time difference of 0.248091 secs
## [1] 78
## Time difference of 0.2544405 secs
## [1] 79
## Time difference of 0.249547 secs
## [1] 80
## Time difference of 0.247648 secs
## [1] 81
## Time difference of 0.2493744 secs
## [1] 82
## Time difference of 0.2482848 secs
## [1] 83
## Time difference of 0.2484884 secs
## [1] 84
## Time difference of 0.2487473 secs
## [1] 85
## Time difference of 0.2500906 secs
## [1] 86
## Time difference of 0.2484577 secs
## [1] 87
## Time difference of 0.2488511 secs
## [1] 88
## Time difference of 0.2483068 secs
## [1] 89
## Time difference of 0.2485266 secs
## [1] 90
## Time difference of 0.249666 secs
## [1] 91
## Time difference of 0.2490311 secs
## [1] 92
## Time difference of 0.2492802 secs
## [1] 93
## Time difference of 0.250612 secs
## [1] 94
## Time difference of 0.2547398 secs
## [1] 95
## Time difference of 0.2492371 secs
## [1] 96
## Time difference of 0.2488213 secs
## [1] 97
## Time difference of 0.2492409 secs
## [1] 98
## Time difference of 0.2482474 secs
## [1] 99
## Time difference of 0.2488101 secs
## [1] 100
## Time difference of 0.2494605 secs
## [1] "Build return matrix"
## [1] "Done"
## Time difference of 25.05278 secs## Formating the data as long ##
long = ensemble_result %>%
as.data.frame() %>%
tibble::rownames_to_column("model") %>%
pivot_longer(2:length(.),names_to = "row_id") %>%
inner_join(gnps_annotations,by = "row_id")
## Extracting top metabolites ##
top = long %>%
group_by(row_id,compound_name_cleaned) %>%
dplyr::summarise(sum = sum(value)) %>%
dplyr::arrange(-sum)
## Pulling top N metabolites ##
topn = top[1:10,] %>%
pull(row_id)
topn_metabolomics = tibble(row_id = topn)
write_csv(topn_metabolomics,"01_r_processed_data/top10_faecalis_v_faecium_metabolomics_biomarkers.csv")
## Filtering the data, combining annotation names ##
filtered = long %>%
filter(row_id %in% topn)
## Plotting the data ##
EFS_plot = filtered %>%
ggplot(aes(reorder(compound_name_cleaned,value),value,fill = model))+
geom_col()+
coord_flip()+
viridis::scale_fill_viridis(discrete = TRUE)+
xlab("Metabolite Biomarkers")+
theme(legend.position = "bottom")
EFS_plot
#ggsave("02_figures/f5_faecalis_v_faecium_metabolomics/f5b.pdf",f5b,height = 8,width = 8)mycophenolic acid and pc(17:1/0:0 + c25h51n1o7p1) were detected as two seperate row IDs. They were both found to be among the best biomarkers independantly from each other.
Plotting
ROC_data = dplyr::filter(normalized_data,row_id %in% topn[1:2]) %>%
dplyr::arrange(factor(row_id,topn)) %>%
inner_join(gnps_annotations,by = "row_id") %>%
filter(condition != "healthy") %>%
mutate(infection_status = condition)
top = make_violin_plot_metabolomics(ROC_data,topn[1:2],ncol = 2,pname = "p.signif")
bottom = make_roc_metabolomics(ROC_data,ncol = 2)
f5b = ggpubr::ggarrange(top,bottom,nrow = 2)
f5b
ggsave("02_figures/f5_faecalis_v_faecium_metabolomics/f5b.pdf",f5b,height = 8, width = 8)Figure 6 - Prediction of Clinical Outcome
A Volcano plot of Proteomics
### Need to rerun msstats, this time changing the condition to be death upon admission ###
# read in clincal metadata
cmd = read_csv("00_r_input_data/clinical_metadata.csv") %>%
select(BioReplicate = sample_id,death_during_admission)
# Read MSstats.csv file.
data = read_csv("00_r_input_data/msstats.csv")
#Note that the condition of the bridge channel needs to be labeled as Norm
annotation = read_csv("00_r_input_data/MSstatsTMT_annotation.csv") %>%
inner_join(cmd,by = "BioReplicate") %>%
select(-Condition) %>%
rename(Condition = death_during_admission) %>%
mutate(Condition = as.factor(as.character(Condition))) %>%
mutate(Condition = case_when(is.na(Condition) ~ "unknown",
TRUE ~ Condition))
mstats = MSstatsTMT::PhilosophertoMSstatsTMTFormat(input = data,annotation)
# use MSstats for protein summarization
quant.msstats = MSstatsTMT::proteinSummarization(mstats,
method="msstats",
global_norm=TRUE,
reference_norm=TRUE,
remove_norm_channel = TRUE,
remove_empty_channel = TRUE)
# Taking the pl_data and formatting appropriately
pl_data = quant.msstats$ProteinLevelData %>%
tidyr::separate(Protein, c("orientation","uniprotid","locus"),sep ="\\|") %>%
# Converting to gene name
inner_join(protein_to_gene,by = c("uniprotid" = "ProteinID")) %>%
mutate(sample_id = paste0(Mixture,".",Channel))
write_csv(pl_data,"01_r_processed_data/quant_msstats_protein_level_data_death.csv")
### Preforming the statistics ###
test.pairwise = MSstatsTMT::groupComparisonTMT(quant.msstats,
moderated = TRUE)
stats =test.pairwise$ComparisonResult %>%
filter(Label == "FALSE vs TRUE") %>%
separate(Protein, c("orientation","uniprotid","locus"),sep ="\\|") %>%
inner_join(protein_to_gene,by = c("uniprotid" = "ProteinID"))
write_csv(stats,"01_r_processed_data/death_vs_life.csv")dvl = read_csv("01_r_processed_data/death_vs_life.csv") %>%
filter(is.na(issue))
f6a = dvl %>%
ggplot(aes(log2FC,-log10(adj.pvalue)))+
geom_point(size = 0.5)+
geom_hline(yintercept = -log10(0.05),linetype = "dashed",color = "red")+
xlab("Log2FC (Live/Dead)")+
ggprism::theme_prism(base_size = 10)
f6a
ggsave("02_figures/f6_prediction_of_clinical_outcome/f6a.pdf",f6a,height = 3.8, width = 4)
# Extracting the signficiant proteins from the volcano plot for later use.
sig = dvl %>%
filter(adj.pvalue <= 0.05)
up_in_live = sig %>%
filter(log2FC>0) %>%
pull(Gene)
up_in_dead = sig %>%
filter(log2FC<0) %>%
pull(Gene)B - GO Enrichment in Dead
#Getting all the go annotations
go = read_csv("00_r_input_data/go_annotations.csv") %>%
dplyr::select(1,length(.)) %>%
tidyr::separate_rows(Gene.Ontology.IDs, sep = ';') %>%
dplyr::mutate(Gene.Ontology.IDs = gsub(" ","",Gene.Ontology.IDs)) %>%
# Going from Protein ID to gene ID
left_join(protein_to_gene,by = c("Entry" = "ProteinID")) %>%
select(-Entry) %>%
select(Entry = Gene,Gene.Ontology.IDs)
#usng enricher
TERM2GENE = go %>%
dplyr::select(TERM = Gene.Ontology.IDs,GENE = Entry) %>%
as.data.frame()
TERM2NAME = AnnotationDbi::select(GO.db::GO.db,keys = unique(TERM2GENE$TERM),columns = c("TERM")) %>%
as.data.frame()
# Perfomring GO enrichment
enrich = clusterProfiler::enricher(up_in_dead,
universe = go$Entry,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME)
f6b = clusterProfiler::cnetplot(enrich, categorySize="pvalue",circular = FALSE,categorySize="pvalue",cex_label_category = 1,showCategory = 100,cex_label_gene = 0.5)+
theme(legend.position = "none")
f6b 
ggsave("02_figures/f6_prediction_of_clinical_outcome/f6b.pdf",f6b,height = 4, width = 4)C - GO Enrichment in Live
#Getting all the go annotations
go = read_csv("00_r_input_data/go_annotations.csv") %>%
dplyr::select(1,length(.)) %>%
tidyr::separate_rows(Gene.Ontology.IDs, sep = ';') %>%
dplyr::mutate(Gene.Ontology.IDs = gsub(" ","",Gene.Ontology.IDs)) %>%
# Going from Protein ID to gene ID
left_join(protein_to_gene,by = c("Entry" = "ProteinID")) %>%
select(-Entry) %>%
select(Entry = Gene,Gene.Ontology.IDs)
#usng enricher
TERM2GENE = go %>%
dplyr::select(TERM = Gene.Ontology.IDs,GENE = Entry) %>%
as.data.frame()
TERM2NAME = AnnotationDbi::select(GO.db::GO.db,keys = unique(TERM2GENE$TERM),columns = c("TERM")) %>%
as.data.frame()
# Perfomring GO enrichment
enrich = clusterProfiler::enricher(up_in_live,
universe = go$Entry,
TERM2GENE = TERM2GENE,
TERM2NAME = TERM2NAME)
f6c = clusterProfiler::cnetplot(enrich, categorySize="pvalue",
circular = FALSE,
categorySize="pvalue",
cex_label_category = 1,
showCategory = 100,
cex_label_gene = 0.5, )+
theme(legend.position = "none")
f6c
ggsave("02_figures/f6_prediction_of_clinical_outcome/f6c.pdf",f6c)D - Proteomic Biomarkers to predict death
Data Processing
##reading in the data##
pl = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv") %>%
rename(tmt_id = sample_id)
## Clinical metadata ##
map = readxl::read_excel("00_r_input_data/proteomics_sample_mapping.xlsx") %>%
mutate(tmt_id = paste0("plex",plex_number,".",TMT_channel)) %>%
select(tmt_id,sample_id = sample_name)
md = read_csv("00_r_input_data/clinical_metadata.csv")
md_final = inner_join(map,md,by = "sample_id")
clean_clinical = inner_join(pl,md_final, by = "tmt_id") %>%
# Don't need the healthy samples here
filter(Condition != "Healthy") %>%
#IGKC gets identified by EFS even though it is a terrible biomarker. Not sure why, just going to remove it
filter(Gene != "IGKC") %>%
# - messes up EFS for some reason, very annoying.
mutate(Gene = gsub("-","_",Gene)) %>%
select(sample_id,death_during_admission,Gene,Abundance) %>%
tidyr::pivot_wider(names_from = Gene,
values_from = Abundance) %>%
#Removing samples with NA.
filter(!is.na(death_during_admission)) %>%
dplyr::select_if(~ !any(is.na(.))) %>%
#needs to be a dataframe or EFS freaks out.
as.data.frame() %>%
select(-sample_id)
set.seed(2)
ensemble_result = EFS::ensemble_fs(clean_clinical,1,cor_threshold = 0)## [1] "default value for NA_threshold = 0.2"
## [1] "default value for runs = 100"
## [1] "default value for selection is c(TRUE, TRUE, TRUE,TRUE, TRUE, TRUE, FALSE, FALSE)"
## [1] "Start Median"
## [1] "Start Pearson"
## [1] "Start Spearman"
## [1] "Start LogReg"
## [1] "Start RF"
## [1] 1
## Time difference of 0.2388246 secs
## [1] 2
## Time difference of 0.2382493 secs
## [1] 3
## Time difference of 0.2379012 secs
## [1] 4
## Time difference of 0.2380173 secs
## [1] 5
## Time difference of 0.2380569 secs
## [1] 6
## Time difference of 0.2392192 secs
## [1] 7
## Time difference of 0.239584 secs
## [1] 8
## Time difference of 0.2403662 secs
## [1] 9
## Time difference of 0.2408121 secs
## [1] 10
## Time difference of 0.2384887 secs
## [1] 11
## Time difference of 0.2390845 secs
## [1] 12
## Time difference of 0.2397776 secs
## [1] 13
## Time difference of 0.2389698 secs
## [1] 14
## Time difference of 0.2411878 secs
## [1] 15
## Time difference of 0.239996 secs
## [1] 16
## Time difference of 0.2404115 secs
## [1] 17
## Time difference of 0.2382421 secs
## [1] 18
## Time difference of 0.2646215 secs
## [1] 19
## Time difference of 0.2389402 secs
## [1] 20
## Time difference of 0.2380698 secs
## [1] 21
## Time difference of 0.2377357 secs
## [1] 22
## Time difference of 0.2391071 secs
## [1] 23
## Time difference of 0.2386279 secs
## [1] 24
## Time difference of 0.2372608 secs
## [1] 25
## Time difference of 0.2371483 secs
## [1] 26
## Time difference of 0.237448 secs
## [1] 27
## Time difference of 0.23721 secs
## [1] 28
## Time difference of 0.2376616 secs
## [1] 29
## Time difference of 0.2381759 secs
## [1] 30
## Time difference of 0.2399197 secs
## [1] 31
## Time difference of 0.2394013 secs
## [1] 32
## Time difference of 0.2391608 secs
## [1] 33
## Time difference of 0.2384644 secs
## [1] 34
## Time difference of 0.2377789 secs
## [1] 35
## Time difference of 0.2380095 secs
## [1] 36
## Time difference of 0.2386534 secs
## [1] 37
## Time difference of 0.2378578 secs
## [1] 38
## Time difference of 0.2374029 secs
## [1] 39
## Time difference of 0.2380579 secs
## [1] 40
## Time difference of 0.2405787 secs
## [1] 41
## Time difference of 0.242275 secs
## [1] 42
## Time difference of 0.2507811 secs
## [1] 43
## Time difference of 0.2398386 secs
## [1] 44
## Time difference of 0.2385688 secs
## [1] 45
## Time difference of 0.2378111 secs
## [1] 46
## Time difference of 0.2384872 secs
## [1] 47
## Time difference of 0.2393737 secs
## [1] 48
## Time difference of 0.2384326 secs
## [1] 49
## Time difference of 0.2368307 secs
## [1] 50
## Time difference of 0.2384396 secs
## [1] 51
## Time difference of 0.2380364 secs
## [1] 52
## Time difference of 0.2378435 secs
## [1] 53
## Time difference of 0.238651 secs
## [1] 54
## Time difference of 0.2375729 secs
## [1] 55
## Time difference of 0.2388303 secs
## [1] 56
## Time difference of 0.2381639 secs
## [1] 57
## Time difference of 0.2383633 secs
## [1] 58
## Time difference of 0.2381697 secs
## [1] 59
## Time difference of 0.2374308 secs
## [1] 60
## Time difference of 0.2383893 secs
## [1] 61
## Time difference of 0.2382855 secs
## [1] 62
## Time difference of 0.2371404 secs
## [1] 63
## Time difference of 0.2378278 secs
## [1] 64
## Time difference of 0.2363172 secs
## [1] 65
## Time difference of 0.2380781 secs
## [1] 66
## Time difference of 0.2396228 secs
## [1] 67
## Time difference of 0.247962 secs
## [1] 68
## Time difference of 0.2391088 secs
## [1] 69
## Time difference of 0.2384813 secs
## [1] 70
## Time difference of 0.2393396 secs
## [1] 71
## Time difference of 0.2400808 secs
## [1] 72
## Time difference of 0.2387977 secs
## [1] 73
## Time difference of 0.2381752 secs
## [1] 74
## Time difference of 0.2384057 secs
## [1] 75
## Time difference of 0.2394662 secs
## [1] 76
## Time difference of 0.2393548 secs
## [1] 77
## Time difference of 0.2375376 secs
## [1] 78
## Time difference of 0.2381115 secs
## [1] 79
## Time difference of 0.2378492 secs
## [1] 80
## Time difference of 0.2371721 secs
## [1] 81
## Time difference of 0.2375021 secs
## [1] 82
## Time difference of 0.2389581 secs
## [1] 83
## Time difference of 0.2378228 secs
## [1] 84
## Time difference of 0.2395983 secs
## [1] 85
## Time difference of 0.2374671 secs
## [1] 86
## Time difference of 0.2368007 secs
## [1] 87
## Time difference of 0.2384114 secs
## [1] 88
## Time difference of 0.2375865 secs
## [1] 89
## Time difference of 0.2377334 secs
## [1] 90
## Time difference of 0.2382908 secs
## [1] 91
## Time difference of 0.2480917 secs
## [1] 92
## Time difference of 0.2388194 secs
## [1] 93
## Time difference of 0.2391233 secs
## [1] 94
## Time difference of 0.2387002 secs
## [1] 95
## Time difference of 0.2379842 secs
## [1] 96
## Time difference of 0.2378576 secs
## [1] 97
## Time difference of 0.2373199 secs
## [1] 98
## Time difference of 0.2387133 secs
## [1] 99
## Time difference of 0.239501 secs
## [1] 100
## Time difference of 0.2372859 secs
## [1] "Build return matrix"
## [1] "Done"
## Time difference of 24.06239 secslong = ensemble_result %>%
as.data.frame() %>%
tibble::rownames_to_column("model") %>%
pivot_longer(2:length(.),names_to = "Gene") %>%
#converting genes back to proper notation.
mutate(Gene = gsub("_","-",Gene))
top = long %>%
group_by(Gene) %>%
dplyr::summarise(sum = sum(value)) %>%
dplyr::arrange(-sum)
topn = top[1:10,] %>%
pull(Gene)
topn_out = tibble(Gene = topn)
write_csv(topn_out,"01_r_processed_data/top10_live_v_dead_proteomics_biomarkers.csv")
filtered = long %>%
filter(Gene %in% topn)
EFS = filtered %>%
ggplot(aes(reorder(Gene,value),value,fill = model))+
geom_col()+
coord_flip()+
viridis::scale_fill_viridis(discrete = TRUE)+
ggtitle("Top 10 Performing Protein Biomarkers")+
xlab("Protein Biomarkers")
EFS
Plotting
source("functions.R")
data = inner_join(pl,md_final, by = "tmt_id") %>%
filter(Condition != "Healthy") %>%
mutate(Condition = death_during_admission) %>%
filter(!is.na(Condition))
top = make_violin_plot_proteomics(data,topn[1:2],pname = "p.signif",ncol = 2)
bottom = make_roc_proteomics(data,topn[1:2],ncol = 2)
f6d = ggpubr::ggarrange(top,bottom,nrow = 2)
f6d
ggsave("02_figures/f6_prediction_of_clinical_outcome/f6d.pdf",f6d,height = 8, width = 8)E- Metabolomics Volcano Plot
## Reading in Data
normalized_data = read_csv("01_r_processed_data/normalized_metabolomics.csv") %>%
filter(!is.na(death_during_admission))
# Performing statistics so we can make a volcano plot
stat = normalized_data %>%
dplyr::group_by(row_id) %>%
rstatix::t_test(norm_value ~ death_during_admission) %>%
rstatix::adjust_pvalue(p.col = "p",output.col = "p.adj_fdr", method = "fdr")
# Calculating log2fc
log2fc_data = normalized_data %>%
dplyr::select(metabolomics_id,row_id,death_during_admission,norm_value)
# Extracting just the faecalis data
dead = log2fc_data %>%
filter(death_during_admission == TRUE) %>%
group_by(row_id) %>%
summarise(dead = mean(norm_value))
# Extracting just the faecium data
live = log2fc_data %>%
filter(death_during_admission == FALSE) %>%
group_by(row_id) %>%
summarise(live = mean(norm_value))
# Calculating log2fc
log2fc = inner_join(dead,live,by = "row_id") %>%
mutate(log2fc_live_dead = log2(live/dead))
# Combining log2fc and statistics
vp = inner_join(stat,log2fc,by = "row_id") %>%
mutate(annotated = case_when(row_id %in% gnps_annotations$row_id ~ TRUE,
TRUE ~ FALSE))
annotations = vp %>%
filter(p.adj_fdr <= 0.05,
annotated == TRUE) %>%
inner_join(gnps_annotations, by = "row_id")
f6e = vp %>%
ggplot(aes(log2fc_live_dead,-log10(p.adj_fdr),color = annotated))+
geom_point(size = 0.5)+
geom_hline(yintercept = -log10(0.05),linetype = "dashed",color = "red")+
xlab("log2FC (Live/Dead)")+
theme(legend.position = "bottom")+
ggrepel::geom_label_repel(data = annotations,aes(label = compound_name_cleaned,x= log2fc_live_dead, y = -log10(p.adj_fdr)),size = 2.5,label.size = 0.1,nudge_y = 0.1)+
ggprism::theme_prism(base_size = 10)+
theme(legend.position = "none")
f6e
ggsave(filename = "02_figures/f6_prediction_of_clinical_outcome/f6e.pdf",f6e,height = 3.8, width = 4)F- Evaluation of Metabolomics Biomarkers to Predict Death
## Reading in the normalized data ##
normalized_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv")
## Switching to the wide format ###
norm_wide = normalized_data %>%
select(-row_m_z,-row_retention_time,-value) %>%
tidyr::pivot_wider(names_from = row_id,values_from = norm_value)
## Looking for only columns that were annotated ##
gnps_and_in_norm = intersect(unique(gnps_annotations$row_id), colnames(norm_wide))
gnps_annotated_metabolites = norm_wide %>%
as.data.frame() %>%
select(death_during_admission,gnps_and_in_norm) %>%
filter(!is.na(death_during_admission))
# Seting seed incase this isn't done internally in EFS
set.seed(2)
## Running EFS feature selection ##
ensemble_result = EFS::ensemble_fs(gnps_annotated_metabolites,
classnumber = 1,cor_threshold = 0)## [1] "default value for NA_threshold = 0.2"
## [1] "default value for runs = 100"
## [1] "default value for selection is c(TRUE, TRUE, TRUE,TRUE, TRUE, TRUE, FALSE, FALSE)"
## [1] "Start Median"
## [1] "Start Pearson"
## [1] "Start Spearman"
## [1] "Start LogReg"
## [1] "Start RF"
## [1] 1
## Time difference of 0.1942923 secs
## [1] 2
## Time difference of 0.1943436 secs
## [1] 3
## Time difference of 0.1948893 secs
## [1] 4
## Time difference of 0.1960897 secs
## [1] 5
## Time difference of 0.1964512 secs
## [1] 6
## Time difference of 0.1960495 secs
## [1] 7
## Time difference of 0.1964924 secs
## [1] 8
## Time difference of 0.1958582 secs
## [1] 9
## Time difference of 0.1948428 secs
## [1] 10
## Time difference of 0.1955628 secs
## [1] 11
## Time difference of 0.1956298 secs
## [1] 12
## Time difference of 0.1962173 secs
## [1] 13
## Time difference of 0.214503 secs
## [1] 14
## Time difference of 0.1947675 secs
## [1] 15
## Time difference of 0.1950436 secs
## [1] 16
## Time difference of 0.193929 secs
## [1] 17
## Time difference of 0.1949337 secs
## [1] 18
## Time difference of 0.1939337 secs
## [1] 19
## Time difference of 0.194304 secs
## [1] 20
## Time difference of 0.1947954 secs
## [1] 21
## Time difference of 0.1949487 secs
## [1] 22
## Time difference of 0.1928654 secs
## [1] 23
## Time difference of 0.1946545 secs
## [1] 24
## Time difference of 0.194489 secs
## [1] 25
## Time difference of 0.1931984 secs
## [1] 26
## Time difference of 0.1944456 secs
## [1] 27
## Time difference of 0.1943586 secs
## [1] 28
## Time difference of 0.1942968 secs
## [1] 29
## Time difference of 0.194309 secs
## [1] 30
## Time difference of 0.1955562 secs
## [1] 31
## Time difference of 0.1965129 secs
## [1] 32
## Time difference of 0.1953988 secs
## [1] 33
## Time difference of 0.1951709 secs
## [1] 34
## Time difference of 0.1947422 secs
## [1] 35
## Time difference of 0.1959519 secs
## [1] 36
## Time difference of 0.1956093 secs
## [1] 37
## Time difference of 0.1956692 secs
## [1] 38
## Time difference of 0.1954944 secs
## [1] 39
## Time difference of 0.1961546 secs
## [1] 40
## Time difference of 0.1958635 secs
## [1] 41
## Time difference of 0.1961136 secs
## [1] 42
## Time difference of 0.2036498 secs
## [1] 43
## Time difference of 0.1949241 secs
## [1] 44
## Time difference of 0.1938779 secs
## [1] 45
## Time difference of 0.1935399 secs
## [1] 46
## Time difference of 0.1941926 secs
## [1] 47
## Time difference of 0.1945207 secs
## [1] 48
## Time difference of 0.1940196 secs
## [1] 49
## Time difference of 0.194773 secs
## [1] 50
## Time difference of 0.1937551 secs
## [1] 51
## Time difference of 0.1936898 secs
## [1] 52
## Time difference of 0.1932113 secs
## [1] 53
## Time difference of 0.1955125 secs
## [1] 54
## Time difference of 0.1944671 secs
## [1] 55
## Time difference of 0.1929574 secs
## [1] 56
## Time difference of 0.1942315 secs
## [1] 57
## Time difference of 0.1933637 secs
## [1] 58
## Time difference of 0.1929708 secs
## [1] 59
## Time difference of 0.1934581 secs
## [1] 60
## Time difference of 0.192693 secs
## [1] 61
## Time difference of 0.1937892 secs
## [1] 62
## Time difference of 0.194901 secs
## [1] 63
## Time difference of 0.1946435 secs
## [1] 64
## Time difference of 0.1944411 secs
## [1] 65
## Time difference of 0.193886 secs
## [1] 66
## Time difference of 0.1946685 secs
## [1] 67
## Time difference of 0.1928484 secs
## [1] 68
## Time difference of 0.1948073 secs
## [1] 69
## Time difference of 0.194073 secs
## [1] 70
## Time difference of 0.1948652 secs
## [1] 71
## Time difference of 0.2044704 secs
## [1] 72
## Time difference of 0.1955822 secs
## [1] 73
## Time difference of 0.1945601 secs
## [1] 74
## Time difference of 0.1935005 secs
## [1] 75
## Time difference of 0.1956577 secs
## [1] 76
## Time difference of 0.1955831 secs
## [1] 77
## Time difference of 0.1944203 secs
## [1] 78
## Time difference of 0.1950929 secs
## [1] 79
## Time difference of 0.1939764 secs
## [1] 80
## Time difference of 0.1941063 secs
## [1] 81
## Time difference of 0.195457 secs
## [1] 82
## Time difference of 0.1935701 secs
## [1] 83
## Time difference of 0.1945205 secs
## [1] 84
## Time difference of 0.1941588 secs
## [1] 85
## Time difference of 0.1937833 secs
## [1] 86
## Time difference of 0.1943376 secs
## [1] 87
## Time difference of 0.1946282 secs
## [1] 88
## Time difference of 0.193825 secs
## [1] 89
## Time difference of 0.1944258 secs
## [1] 90
## Time difference of 0.1951914 secs
## [1] 91
## Time difference of 0.1946559 secs
## [1] 92
## Time difference of 0.1942747 secs
## [1] 93
## Time difference of 0.1943419 secs
## [1] 94
## Time difference of 0.19487 secs
## [1] 95
## Time difference of 0.1941037 secs
## [1] 96
## Time difference of 0.192903 secs
## [1] 97
## Time difference of 0.1946273 secs
## [1] 98
## Time difference of 0.1949513 secs
## [1] 99
## Time difference of 0.1952524 secs
## [1] 100
## Time difference of 0.2035351 secs
## [1] "Build return matrix"
## [1] "Done"
## Time difference of 19.61908 secs## Formating the data as long ##
long = ensemble_result %>%
as.data.frame() %>%
tibble::rownames_to_column("model") %>%
pivot_longer(2:length(.),names_to = "row_id")
## Extracting top metabolites ##
top = long %>%
group_by(row_id) %>%
dplyr::summarise(sum = sum(value)) %>%
dplyr::arrange(-sum)
## Pulling top N metabolites ##
topn = top[1:10,] %>%
pull(row_id)
topn_out = tibble(row_id = topn )
write_csv(topn_out,"01_r_processed_data/top10_live_v_dead_metabolomics_biomarkers.csv")
data = normalized_data %>%
inner_join(gnps_annotations,by = "row_id") %>%
filter(!is.na(death_during_admission)) %>%
filter(row_id %in% topn[1:2]) %>%
mutate(condition = as.character(death_during_admission))
top = make_violin_plot_metabolomics(data,ncol = 2,topn = topn[1:2],pname = "p.signif")
bottom = make_roc_metabolomics(data,ncol = 2)
f6f = ggpubr::ggarrange(top,bottom,nrow = 2)
f6f
ggsave("02_figures/f6_prediction_of_clinical_outcome/f6f.pdf",f6f,height = 8,width = 8)Suplemental Figures
S1 - Overview of Analysis Workflow
This is generated in draw.io
S2 - Networks
See Figure 2 B and Figure 3 C
S3 - Known Biomarkers of Infection
A - ROC curves.
pl = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
known_indicatiors_of_bacterial_infection = c("CRP","SAA1")
f = dplyr::filter(pl,Gene %in% known_indicatiors_of_bacterial_infection) %>%
mutate(Condition = case_when(Condition == "Healthy" ~ "Healthy",
TRUE ~ "Infected"))
data = f
list = list()
for(i in known_indicatiors_of_bacterial_infection){
loop = filter(data, Gene == i)
mod = glm(as.factor(Condition) ~ Abundance,
data = data,
family = "binomial")
predicted <- predict(mod ,family = "binomial", loop, type="response")
#define object to plot
rocobj <- pROC::roc(loop$Condition, predicted)
auc = round(as.numeric(gsub(".*: ", "",rocobj$auc)),2)
#create ROC plot
p1 = pROC::ggroc(rocobj)+
ggtitle(paste0(i," AUC: ",auc))
list[[i]] = p1
}
s3a = do.call("grid.arrange", c(list, ncol = 2),)
s3a## TableGrob (1 x 2) "arrange": 2 grobs
## z cells name grob
## CRP 1 (1-1,1-1) arrange gtable[layout]
## SAA1 2 (1-1,2-2) arrange gtable[layout]ggsave("03_supplemental/s3_known_biomarkers_of_infection/s3a.pdf",s3a)B - Violin Plot
#Loading and processing the data
#reading in the data
pl = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
# C-reactive Protein = "P02741
# Calcitonin = P01258
# Serum Amyloid A = P0DJI8
known_indicatiors_of_bacterial_infection = c("CRP","SAA1")
f = dplyr::filter(pl,Gene %in% known_indicatiors_of_bacterial_infection)
library(rstatix)
library(ggpubr)
ttest <-
f %>%
group_by(Gene) %>%
rstatix::t_test(Abundance ~ Condition) %>%
add_y_position()
s3b = f %>%
ggplot(aes(Condition,Abundance))+
geom_violin()+
geom_point()+
ggpubr::stat_pvalue_manual(ttest)+
ylab("Abundance")+
facet_wrap(~Gene)
s3b
ggsave("03_supplemental/s3_known_biomarkers_of_infection/s3b.pdf",s3b)S4 - Top 10 identified biomarkers Healthy vs Infected
A- proteomics violin plot
# Reading in data
protein_sum = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
topn = read_csv("01_r_processed_data/top10_infected_v_healthy_proteomics_biomarkers.csv") %>% pull(Gene)
# Processing data
ROC_data = dplyr::filter(protein_sum,Gene %in% topn) #%>%
#dplyr::mutate(Condition = case_when(Condition %in% c("Faecalis","Faecium") ~ "infected",
# TRUE ~ Condition))
source("functions.R")
s4a = make_violin_plot_proteomics(ROC_data,topn)
s4a## TableGrob (2 x 5) "arrange": 10 grobs
## z cells name grob
## SERPINA3 1 (1-1,1-1) arrange gtable[layout]
## LRG1 2 (1-1,2-2) arrange gtable[layout]
## CLEC3B 3 (1-1,3-3) arrange gtable[layout]
## TF 4 (1-1,4-4) arrange gtable[layout]
## APOA2 5 (1-1,5-5) arrange gtable[layout]
## LBP 6 (2-2,1-1) arrange gtable[layout]
## ENO1 7 (2-2,2-2) arrange gtable[layout]
## GSN 8 (2-2,3-3) arrange gtable[layout]
## SERPINA5 9 (2-2,4-4) arrange gtable[layout]
## FGA 10 (2-2,5-5) arrange gtable[layout]ggsave("03_supplemental/s4_biomarkers_healthy_v_infected/s4a.pdf",s4a, height = 6 ,width=16)B- Proteomics ROC
# Processing data
ROC_data = dplyr::filter(protein_sum,Gene %in% topn) %>%
dplyr::mutate(Condition = case_when(Condition %in% c("Faecalis","Faecium") ~ "infected",
TRUE ~ Condition))
source("functions.R")
s4b = make_roc_proteomics(ROC_data,topn)
s4b## TableGrob (2 x 5) "arrange": 10 grobs
## z cells name grob
## SERPINA3 1 (1-1,1-1) arrange gtable[layout]
## LRG1 2 (1-1,2-2) arrange gtable[layout]
## CLEC3B 3 (1-1,3-3) arrange gtable[layout]
## TF 4 (1-1,4-4) arrange gtable[layout]
## APOA2 5 (1-1,5-5) arrange gtable[layout]
## LBP 6 (2-2,1-1) arrange gtable[layout]
## ENO1 7 (2-2,2-2) arrange gtable[layout]
## GSN 8 (2-2,3-3) arrange gtable[layout]
## SERPINA5 9 (2-2,4-4) arrange gtable[layout]
## FGA 10 (2-2,5-5) arrange gtable[layout]ggsave("03_supplemental/s4_biomarkers_healthy_v_infected/s4b.pdf",s4b, height = 6 ,width=16)C - Metabolomics Violin
## Reading in the normalized data ##
normalized_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv")
topn = read_csv("01_r_processed_data/top10_infected_v_healthy_metabolomics_biomarkers.csv") %>% pull(row_id)
ROC_data = dplyr::filter(normalized_data,row_id %in% topn) %>%
dplyr::arrange(factor(row_id,topn)) %>%
inner_join(gnps_annotations,by = c("row_id" = "row_id"))
source("functions.R")
s4c = make_violin_plot_metabolomics(ROC_data,topn)
ggsave("03_supplemental/s4_biomarkers_healthy_v_infected/s4c.pdf",s4c,height = 6 ,width=16)D - Metabolomics ROC
## Reading in the normalized data ##
normalized_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv")
topn = read_csv("01_r_processed_data/top10_infected_v_healthy_metabolomics_biomarkers.csv") %>% pull(row_id)
ROC_data = dplyr::filter(normalized_data,row_id %in% topn) %>%
dplyr::arrange(factor(row_id,topn)) %>%
inner_join(gnps_annotations,by = c("row_id" = "row_id")) %>%
dplyr::mutate(condition = case_when(condition %in% c("faecalis","faecium") ~ "infected",
TRUE ~ condition))
source("functions.R")
s4d = make_roc_metabolomics(ROC_data,)
ggsave("03_supplemental/s4_biomarkers_healthy_v_infected/s4d.pdf",s4d,height = 6 ,width=16)S5 - Top 10 Faecalis vs Faecium
A Proteomics violin
# Reading in data
protein_sum = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
topn = read_csv("01_r_processed_data/top10_faecalis_v_faecium_proteomics_biomarkers.csv") %>% pull(Gene)
# Processing data
ROC_data = dplyr::filter(protein_sum,Gene %in% topn)
source("functions.R")
s5a = make_violin_plot_proteomics(ROC_data,topn)
s5a## TableGrob (2 x 5) "arrange": 10 grobs
## z cells name grob
## IGKV2-30 1 (1-1,1-1) arrange gtable[layout]
## RBP4 2 (1-1,2-2) arrange gtable[layout]
## APOC3 3 (1-1,3-3) arrange gtable[layout]
## IGHV6-1 4 (1-1,4-4) arrange gtable[layout]
## CLU 5 (1-1,5-5) arrange gtable[layout]
## PROC 6 (2-2,1-1) arrange gtable[layout]
## SERPINC1 7 (2-2,2-2) arrange gtable[layout]
## IGLV1-47 8 (2-2,3-3) arrange gtable[layout]
## PCYOX1 9 (2-2,4-4) arrange gtable[layout]
## CPN1 10 (2-2,5-5) arrange gtable[layout]ggsave("03_supplemental/s5_biomarkers_faecalis_v_faecium/s5a.pdf",s5a,height = 6 ,width=16)B- Proteomics ROC
# Processing data
ROC_data = dplyr::filter(protein_sum,Gene %in% topn) %>%
filter(Condition != "Healthy")
source("functions.R")
s5b = make_roc_proteomics(ROC_data,topn)
s5b## TableGrob (2 x 5) "arrange": 10 grobs
## z cells name grob
## IGKV2-30 1 (1-1,1-1) arrange gtable[layout]
## RBP4 2 (1-1,2-2) arrange gtable[layout]
## APOC3 3 (1-1,3-3) arrange gtable[layout]
## IGHV6-1 4 (1-1,4-4) arrange gtable[layout]
## CLU 5 (1-1,5-5) arrange gtable[layout]
## PROC 6 (2-2,1-1) arrange gtable[layout]
## SERPINC1 7 (2-2,2-2) arrange gtable[layout]
## IGLV1-47 8 (2-2,3-3) arrange gtable[layout]
## PCYOX1 9 (2-2,4-4) arrange gtable[layout]
## CPN1 10 (2-2,5-5) arrange gtable[layout]ggsave("03_supplemental/s5_biomarkers_faecalis_v_faecium/s5b.pdf",s5b,height = 6 ,width=16)C- metabolomics violin
## Reading in the normalized data ##
normalized_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv")
topn = read_csv("01_r_processed_data/top10_faecalis_v_faecium_metabolomics_biomarkers.csv") %>% pull(row_id)
ROC_data = dplyr::filter(normalized_data,row_id %in% topn) %>%
dplyr::arrange(factor(row_id,topn)) %>%
inner_join(gnps_annotations,by = c("row_id" = "row_id"))
s5c = make_violin_plot_metabolomics(ROC_data,topn)
ggsave("03_supplemental/s5_biomarkers_faecalis_v_faecium/s5c.pdf",s5c,height = 6 ,width=16)D- metabolomics ROC
ROC_data = dplyr::filter(normalized_data,row_id %in% topn) %>%
dplyr::arrange(factor(row_id,topn)) %>%
inner_join(gnps_annotations,by = c("row_id" = "row_id")) %>%
filter(condition != "healthy")
s5d = make_roc_metabolomics(ROC_data)
ggsave("03_supplemental/s5_biomarkers_faecalis_v_faecium/s5d.pdf",s5d,height = 6 ,width=16)S6- Top 10 Mortality vs Survival
A- Proteomics violin
##reading in the data##
pl = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv") %>%
rename(tmt_id = sample_id)
## Clinical metadata ##
map = readxl::read_excel("00_r_input_data/proteomics_sample_mapping.xlsx") %>%
mutate(tmt_id = paste0("plex",plex_number,".",TMT_channel)) %>%
select(tmt_id,sample_id = sample_name)
md = read_csv("00_r_input_data/clinical_metadata.csv")
md_final = inner_join(map,md,by = "sample_id")
clean_clinical = inner_join(pl,md_final, by = "tmt_id") %>%
# Don't need the healthy samples here
filter(Condition != "Healthy") %>%
# - messes up EFS for some reason, very annoying.
mutate(Gene = gsub("-","_",Gene)) %>%
select(sample_id,death_during_admission,Gene,Abundance) %>%
tidyr::pivot_wider(names_from = Gene,
values_from = Abundance) %>%
#Removing samples with NA.
filter(!is.na(death_during_admission)) %>%
dplyr::select_if(~ !any(is.na(.))) %>%
#needs to be a dataframe or EFS freaks out.
as.data.frame() %>%
select(-sample_id)
topn = read_csv("01_r_processed_data/top10_live_v_dead_proteomics_biomarkers.csv") %>%
pull(Gene)
data = inner_join(pl,md_final, by = "tmt_id") %>%
filter(Condition != "Healthy") %>%
mutate(Condition = death_during_admission) %>%
filter(!is.na(Condition))
s6a = make_violin_plot_proteomics(data,topn,pname="p.signif")
ggsave("03_supplemental/s6_biomarkers_survival_v_mortality/s6a.pdf",s6a,height = 6, width = 16)B- Proteomics ROC
s6b = make_roc_proteomics(data,topn)
ggsave("03_supplemental/s6_biomarkers_survival_v_mortality/s6b.pdf",s6b,height = 6, width = 16)C - Metabolomics violin
## Reading in the normalized data ##
normalized_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv")
topn = read_csv("01_r_processed_data/top10_live_v_dead_metabolomics_biomarkers.csv") %>%
pull(row_id)
data = normalized_data %>%
inner_join(gnps_annotations,by = "row_id") %>%
filter(!is.na(death_during_admission)) %>%
filter(row_id %in% topn) %>%
mutate(condition = as.character(death_during_admission))
s6c = make_violin_plot_metabolomics(data,topn,ncol = 5,pname="p.signif")
ggsave("03_supplemental/s6_biomarkers_survival_v_mortality/s6c.pdf",s6c,height = 6,width = 16)D - Metabolomics ROC
s6d = make_roc_metabolomics(data)
ggsave("03_supplemental/s6_biomarkers_survival_v_mortality/s6d.pdf",s6d,height = 6,width = 16)S7- metadata assesment of all biomarkers ided in the study
A- Faecalis vs Faecium Proteomics.
clinical_md = read_csv("00_r_input_data/clinical_metadata.csv")
protein_sum = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
# Defining the input metadata. Need to chose which columns to keep here.
md_assesment_data= protein_sum %>%
filter(Condition != "Healthy") %>%
inner_join(clinical_md, by = c("BioReplicate" = "sample_id")) %>%
select(sample_id,Gene,bmi,bacteremia_duration,charleson_index,congestive_heart_failure,
death_during_admission,history_of_myocardial_infarction,transplant_type,pathogen,icu_admission,
hypotension,gender,race,smoking,polymicrobial_bacteremia,neutropenia_anc_1500,leukemia_or_lymphoma,Abundance) %>%
ungroup()
topn = read_csv("01_r_processed_data/top10_faecalis_v_faecium_proteomics_biomarkers.csv") %>%
pull(Gene)
source("functions.R")
stat = metadata_assesment(md_assesment_data,
metric = Gene,
sample_id = sample_id,
features = topn,
outcome = Abundance,
outcome_c = "Abundance")
stat_final = rstatix::adjust_pvalue(stat,"p",output.col = "p.adj",method = "fdr")
## Plotting ###
s7a = stat_final %>%
ggplot(aes(reorder(variable,p.adj),-log10(p.adj)))+
geom_col()+
geom_hline(yintercept = -log10(0.05),linetype = "dashed", color = "red")+
coord_flip()+
facet_wrap(~Gene,nrow = 2) +
xlab("metadata_field")
s7a
ggsave("03_supplemental/s7_biomarkers_clinical_metadata_assessment/s7a.pdf",s7a,height = 6, width = 16)B- Faecalis vs Faecium Metabolomics.
metabolomics_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv")
m = inner_join(metabolomics_data,gnps_annotations,by = "row_id")
# Defining the input metadata. Need to chose which columns to keep here.
md_assesment_data= m %>%
filter(condition != "healthy") %>%
select(sample_id,row_id,bmi,bacteremia_duration,charleson_index,congestive_heart_failure,
death_during_admission,history_of_myocardial_infarction,transplant_type,pathogen,icu_admission,
hypotension,gender,race,smoking,polymicrobial_bacteremia,neutropenia_anc_1500,leukemia_or_lymphoma,norm_value) %>%
ungroup()
topn = read_csv("01_r_processed_data/top10_faecalis_v_faecium_metabolomics_biomarkers.csv") %>%
pull(row_id)
source("functions.R")
stat = metadata_assesment(md_assesment_data,
metric = row_id,
sample_id = sample_id,
features = topn,
outcome = norm_value,
outcome_c = "norm_value")
stat_final = rstatix::adjust_pvalue(stat,"p",output.col = "p.adj",method = "fdr")
## Plotting ###
s7b = stat_final %>%
inner_join(gnps_annotations,by = "row_id") %>%
ggplot(aes(reorder(variable,p.adj),-log10(p.adj)))+
geom_col()+
geom_hline(yintercept = -log10(0.05),linetype = "dashed", color = "red")+
coord_flip()+
facet_wrap(~compound_name_cleaned,nrow = 2) +
xlab("metadata_field")
s7b
ggsave("03_supplemental/s7_biomarkers_clinical_metadata_assessment//s7b.pdf",s7b,height = 6, width = 16)C - Proteomic Biomarkers Identified in Live vs Dead.
clinical_md = read_csv("00_r_input_data/clinical_metadata.csv")
protein_sum = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
# Defining the input metadata. Need to chose which columns to keep here.
md_assesment_data= protein_sum %>%
inner_join(clinical_md, by = c("BioReplicate" = "sample_id")) %>%
filter(!is.na(death_during_admission)) %>%
select(sample_id,Gene,bmi,bacteremia_duration,charleson_index,congestive_heart_failure,
death_during_admission,history_of_myocardial_infarction,transplant_type,pathogen,icu_admission,
hypotension,gender,race,smoking,polymicrobial_bacteremia,neutropenia_anc_1500,leukemia_or_lymphoma,Abundance) %>%
ungroup()
topn = read_csv("01_r_processed_data/top10_live_v_dead_proteomics_biomarkers.csv") %>%
pull(Gene)
source("functions.R")
stat = metadata_assesment(md_assesment_data,
metric = Gene,
sample_id = sample_id,
features = topn,
outcome = Abundance,
outcome_c = "Abundance")
stat_final = rstatix::adjust_pvalue(stat,"p",output.col = "p.adj",method = "fdr")
## Plotting ###
s7c = stat_final %>%
ggplot(aes(reorder(variable,p.adj),-log10(p.adj)))+
geom_col()+
geom_hline(yintercept = -log10(0.05),linetype = "dashed", color = "red")+
coord_flip()+
facet_wrap(~Gene,nrow = 2) +
xlab("metadata_field")
s7c
ggsave("03_supplemental/s7_biomarkers_clinical_metadata_assessment//s7c.pdf",s7c,height = 6, width = 16)D - Metabolomic Biomarker Identified in Live vs Dead
metabolomics_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv")
m = inner_join(metabolomics_data,gnps_annotations,by = "row_id")
# Defining the input metadata. Need to chose which columns to keep here.
md_assesment_data= m %>%
filter(!is.na(death_during_admission)) %>%
select(sample_id,row_id,bmi,bacteremia_duration,charleson_index,congestive_heart_failure,
death_during_admission,history_of_myocardial_infarction,transplant_type,pathogen,icu_admission,
hypotension,gender,race,smoking,polymicrobial_bacteremia,neutropenia_anc_1500,leukemia_or_lymphoma,norm_value) %>%
ungroup()
topn = read_csv("01_r_processed_data/top10_live_v_dead_metabolomics_biomarkers.csv") %>%
pull(row_id)
source("functions.R")
stat = metadata_assesment(md_assesment_data,
metric = row_id,
sample_id = sample_id,
features = topn,
outcome = norm_value,
outcome_c = "norm_value")
stat_final = rstatix::adjust_pvalue(stat,"p",output.col = "p.adj",method = "fdr")
## Plotting ###
s7d = stat_final %>%
inner_join(gnps_annotations,by = "row_id") %>%
ggplot(aes(reorder(variable,p.adj),-log10(p.adj)))+
geom_col()+
geom_hline(yintercept = -log10(0.05),linetype = "dashed", color = "red")+
coord_flip()+
facet_wrap(~compound_name_cleaned,nrow = 2) +
xlab("metadata_field")
s7d
ggsave("03_supplemental/s7_biomarkers_clinical_metadata_assessment/s7d.pdf",s7d,height = 6, width = 16)S8 -Cytokine Inference.
Network-based Cytokine Inference Knowledge-based networks were used to infer the relative contributions of major cytokines on the observed proteomic alterations. First, a list of cytokines was manually curated (TGFb, TNF, IFN, IL1-40, CXCL1-16, CCL1-27) and submitted with the significantly altered proteins (uncorrected ANOVA p < 0.05) to the STRING-db tool (all active interaction sources, interaction score > 0.4). The network was exported into a simple tabular output and cytokines were profiled for their enrichment in mortality clusters (pMortality+, pMortality and pMortality) using the following approach. First, cytokines were filtered to have at least five connections to mortality networks. Then, the proportion of each cytokine within the mortality clusters (ie. number of cytokine-mortality cluster connections relative to the total connections within mortality clusters) was compared to the proportion of each cytokine in the entire dataset (ie. number of cytokine-protein connections relative to total connections within the entire dataset). Significance was determined using Chi-square tests and a fold-change was calculated by dividing the proportion of cytokine-mortality connections by the proportion of cytokine connections within the total dataset. A final enrichment score was calculated by multiplying the -log10(p value) from the Chisquare test by the log2(fold-change of enrichment). By including the total number of connections that each cytokine had in the entire dataset as a background, any enrichment bias due to higher annotation rates for popular cytokines in STRING-db is controlled for.
Upstream Data Processing
### Loading data ###
faecalis_v_faecium = readr::read_csv("01_r_processed_data/faecalis_v_faecium.csv")
# Name Mapping
annotations = read_csv("00_r_input_data/gene_uniprot_annotations.csv")
#changing the names to uniprotid
faecalis_v_faecium = left_join(faecalis_v_faecium,annotations,by = c("uniprotid"))
# reading in the curated cytokine list
cytokine_list = read_csv("00_r_input_data/jacobs_cytokine_list.csv",skip = 1)
### Cytokine Network generated for all proteins detected in the study ###
background = faecalis_v_faecium %>%
pull(gene_name)
#exporting file that we can use with the stringdb web tool
#https://www.string-db.org/cgi/input?sessionId=bjuiDWoZnTQb&input_page_show_search=on
all_proteins_and_cytokines = c(background,cytokine_list$cytokines)
write_lines(all_proteins_and_cytokines,"01_r_processed_data/cytokine_inference/00_all_proteins_and_cytokines.csv")
# Results were exported from the stringdb web tool.
# Settings were set as follows: minimum required interaction score = 0.400 (medium), active interaction sources = all (textmining, neighborhood, experiments, gene fusions, databases, co-occurenece, co-expression)
# Gene names that did not find a map were disregarded.
# Was exported as tabular text output with reciprocal edges
# Reading in results input into string db using faecalis_out proteins
string_db_results = read_tsv("01_r_processed_data/cytokine_inference/01_all_proteins_and_cytokines_string_db_out.tsv") %>%
#only keeping connections that are between a protein identified and cytokines
filter(`#node1` %in% background & `node2` %in% cytokine_list$cytokines)
write_csv(string_db_results,"01_r_processed_data/cytokine_inference/02_all_proteins_and_cytokines_string_db_out_filtered.csv")
# counting the number of times a cytokine interacted with a detected protein
n_background_proteins = length(background)
background_cytokines = string_db_results %>%
group_by(node2) %>%
summarise(n = n()) %>%
mutate(condition = "background_present")
write_lines(cytokine_list$cytokines,"01_r_processed_data/cytokine_inference/cytokine_node_list")A - Overview of entire network
This figure is generated in cytoscape using 02_all_proteins_and_cytokines_string_db_out_filtered.csv as the input network. Biopax style was used. Labels were removed. Background was changed to white and the edges were changed to black with a transparency of 75. Nodes in node 2 were selected using the cytokine_node_list file and their size (40), color (dark red) and the shape (diamond) was changed to distinguish the cytokines from other nodes.
4 cytoscape diagrams were created using the initial plot as a template in order to highlight the protein nodes by the protein belonging to each condition. Keeping with the scheme of the paper, healthy samples were colored yellow, faecalis infected samples were colored navy, and faecium colored samples were colored turquoise.
B - Explanation of the stats used to identify signficant cytokines
This figure was created in biorender.
C - Cytokines that were found to be significantly altered among the conditions studied.
Infected
## Doing this for Infected vs Background ##
### Now lets try doing this in Infected vs Healthy ###
infected_vs_healthy = read_csv("01_r_processed_data/Infected_vs_Healthy.csv")
# Name Mapping
annotations = read_csv("00_r_input_data/gene_uniprot_annotations.csv")
#changing the names to uniprotid
infected_vs_healthy = left_join(infected_vs_healthy,annotations,by = c("uniprotid"))
# figuring otu what genes were up in infected
up_in_infected = infected_vs_healthy %>%
filter(log2FC>0) %>%
filter(adj.pvalue <= 0.05) %>%
arrange(adj.pvalue) %>%
pull(gene_name)
n_infected_proteins = length(up_in_infected)
infected_cytokines = string_db_results %>%
filter(`#node1` %in% up_in_infected) %>%
group_by(node2) %>%
summarise(n= n()) %>%
mutate(condition = "infected_present")
background_infected = bind_rows(background_cytokines,infected_cytokines) %>%
pivot_wider(names_from = condition, values_from = n) %>%
mutate(infected_present = replace_na(infected_present,0),
background_absent = n_background_proteins- background_present,
infected_absent = n_infected_proteins - infected_present)
# Making a contingency table
out = data.frame()
for(i in unique(background_infected$node2)){
contingency_table_present = background_infected %>%
filter(node2 == i) %>%
select(infected_present,background_present) %>%
t() %>%
as.data.frame() %>%
rename(present = V1)
contingency_table_absent =background_infected %>%
filter(node2 == i) %>%
select(infected_absent,background_absent) %>%
t() %>%
as.data.frame() %>%
rename(absent = V1)
contingency_table = bind_cols(contingency_table_present,contingency_table_absent) %>%
as.matrix()
stat = rstatix::fisher_test(contingency_table) %>%
mutate(cytokine = i )
out = rbind(stat,out)
}
corrected_pvalues = out %>%
rstatix::adjust_pvalue(method = "fdr")
# Now calculating fold change
log2fc = background_infected %>%
select(node2,
background_present,
infected_present) %>%
mutate(n_background_proteins = n_background_proteins,
n_infected_proteins = n_infected_proteins,
proportion_infected = infected_present/n_infected_proteins,
proportion_background = background_present/n_background_proteins,
log2fc = log2(proportion_infected/proportion_background)) %>%
select(node2,log2fc,proportion_in_condition = proportion_infected,
proportion_in_background = proportion_background)
# all metrics
infected_metrics = corrected_pvalues %>%
inner_join(log2fc,by = c("cytokine" ="node2")) %>%
mutate(enrichment_score = -log10(p.adj) * log2fc) %>%
mutate(condition = "infected")
write_lines(up_in_infected,"01_r_processed_data/cytokine_inference/infected_protein_list")Healthy
### Now lets try doing this in Infected vs Healthy ###
infected_vs_healthy = read_csv("01_r_processed_data/Infected_vs_Healthy.csv")
# Name Mapping
annotations = read_csv("00_r_input_data/gene_uniprot_annotations.csv")
#changing the names to uniprotid
infected_vs_healthy = left_join(infected_vs_healthy,annotations,by = c("uniprotid"))
# up in healthy
up_in_healthy = infected_vs_healthy %>%
filter(log2FC<0) %>%
filter(adj.pvalue <= 0.05) %>%
arrange(adj.pvalue) %>%
pull(gene_name)
n_healthy_proteins = length(up_in_healthy)
healthy_cytokines = string_db_results %>%
filter(`#node1` %in% up_in_healthy) %>%
group_by(node2) %>%
summarise(n= n()) %>%
mutate(condition = "healthy_present")
background_healthy = bind_rows(background_cytokines,healthy_cytokines) %>%
pivot_wider(names_from = condition, values_from = n) %>%
mutate(healthy_present = replace_na(healthy_present,0),
background_absent = n_background_proteins- background_present,
healthy_absent = n_healthy_proteins - healthy_present)
# Making a contingency table
out = data.frame()
for(i in unique(background_healthy$node2)){
contingency_table_present = background_healthy %>%
filter(node2 == i) %>%
select(healthy_present,background_present) %>%
t() %>%
as.data.frame() %>%
rename(present = V1)
contingency_table_absent =background_healthy %>%
filter(node2 == i) %>%
select(healthy_absent,background_absent) %>%
t() %>%
as.data.frame() %>%
rename(absent = V1)
contingency_table = bind_cols(contingency_table_present,contingency_table_absent) %>%
as.matrix()
stat = rstatix::fisher_test(contingency_table) %>%
mutate(cytokine = i )
out = rbind(stat,out)
}
corrected_pvalues = out %>%
rstatix::adjust_pvalue(method = "fdr")
# Now calculating fold change
log2fc = background_healthy %>%
select(node2,
background_present,
healthy_present) %>%
mutate(n_background_proteins = n_background_proteins,
n_healthy_proteins = n_healthy_proteins,
proportion_healthy = healthy_present/n_healthy_proteins,
proportion_background = background_present/n_background_proteins,
log2fc = log2(proportion_healthy/proportion_background)) %>%
select(node2,log2fc,proportion_in_condition = proportion_healthy,
proportion_in_background = proportion_background)
# all metrics
healthy_metrics = corrected_pvalues %>%
inner_join(log2fc,by = c("cytokine" ="node2")) %>%
mutate(enrichment_score = -log10(p.adj) * log2fc,
condition = "healthy")
write_lines(up_in_healthy,"01_r_processed_data/cytokine_inference/healthy_protein_list")Faecalis
### Now lets try doing this in Facalis vs Faecium ###
faecalis_v_faecium = readr::read_csv("01_r_processed_data/faecalis_v_faecium.csv")
# Name Mapping
annotations = read_csv("00_r_input_data/gene_uniprot_annotations.csv")
#changing the names to uniprotid
faecalis_v_faecium = left_join(faecalis_v_faecium,annotations,by = c("uniprotid"))
# up in faecalis
up_in_faecalis = faecalis_v_faecium %>%
filter(log2FC>0) %>%
filter(adj.pvalue <= 0.05) %>%
arrange(adj.pvalue) %>%
pull(gene_name)
n_faecalis_proteins = length(up_in_faecalis)
faecalis_cytokines = string_db_results %>%
filter(`#node1` %in% up_in_faecalis) %>%
group_by(node2) %>%
summarise(n= n()) %>%
mutate(condition = "faecalis_present")
background_faecalis = bind_rows(background_cytokines,faecalis_cytokines) %>%
pivot_wider(names_from = condition, values_from = n) %>%
mutate(faecalis_present = replace_na(faecalis_present,0),
background_absent = n_background_proteins- background_present,
faecalis_absent = n_faecalis_proteins - faecalis_present)
# Making a contingency table
out = data.frame()
for(i in unique(background_faecalis$node2)){
contingency_table_present = background_faecalis %>%
filter(node2 == i) %>%
select(faecalis_present,background_present) %>%
t() %>%
as.data.frame() %>%
rename(present = V1)
contingency_table_absent = background_faecalis %>%
filter(node2 == i) %>%
select(faecalis_absent,background_absent) %>%
t() %>%
as.data.frame() %>%
rename(absent = V1)
contingency_table = bind_cols(contingency_table_present,contingency_table_absent) %>%
as.matrix()
stat = rstatix::fisher_test(contingency_table) %>%
mutate(cytokine = i )
out = rbind(stat,out)
}
corrected_pvalues = out %>%
rstatix::adjust_pvalue(method = "fdr")
# Now calculating fold change
log2fc = background_faecalis %>%
select(node2,
background_present,
faecalis_present) %>%
mutate(n_background_proteins = n_background_proteins,
n_faecalis_proteins = n_faecalis_proteins,
proportion_faecalis = faecalis_present/n_faecalis_proteins,
proportion_background = background_present/n_background_proteins,
log2fc = log2(proportion_faecalis/proportion_background)) %>%
select(node2,log2fc,proportion_in_condition = proportion_faecalis,
proportion_in_background = proportion_background)
# all metrics
faecalis_metrics = corrected_pvalues %>%
inner_join(log2fc,by = c("cytokine" ="node2")) %>%
mutate(enrichment_score = -log10(p.adj) * log2fc,
condition = "faecalis")
write_lines(up_in_faecalis,"01_r_processed_data/cytokine_inference/faecalis_protein_list")Faecium
### Now lets try doing this in Facalis vs Faecium ###
faecalis_v_faecium = readr::read_csv("01_r_processed_data/faecalis_v_faecium.csv")
# Name Mapping
annotations = read_csv("00_r_input_data/gene_uniprot_annotations.csv")
#changing the names to uniprotid
faecalis_v_faecium = left_join(faecalis_v_faecium,annotations,by = c("uniprotid"))
# up in faceium
up_in_faecium = faecalis_v_faecium %>%
filter(log2FC<0) %>%
filter(adj.pvalue <= 0.05) %>%
arrange(adj.pvalue) %>%
pull(gene_name)
n_faecium_proteins = length(up_in_faecium)
faecium_cytokines = string_db_results %>%
filter(`#node1` %in% up_in_faecium) %>%
group_by(node2) %>%
summarise(n= n()) %>%
mutate(condition = "faecium_present")
background_faecium = bind_rows(background_cytokines,faecium_cytokines) %>%
pivot_wider(names_from = condition, values_from = n) %>%
mutate(faecium_present = replace_na(faecium_present,0),
background_absent = n_background_proteins- background_present,
faecium_absent = n_faecium_proteins - faecium_present)
# Making a contingency table
out = data.frame()
for(i in unique(background_faecalis$node2)){
contingency_table_present = background_faecium %>%
filter(node2 == i) %>%
select(faecium_present,background_present) %>%
t() %>%
as.data.frame() %>%
rename(present = V1)
contingency_table_absent = background_faecium %>%
filter(node2 == i) %>%
select(faecium_absent,background_absent) %>%
t() %>%
as.data.frame() %>%
rename(absent = V1)
contingency_table = bind_cols(contingency_table_present,contingency_table_absent) %>%
as.matrix()
stat = rstatix::fisher_test(contingency_table) %>%
mutate(cytokine = i )
out = rbind(stat,out)
}
corrected_pvalues = out %>%
rstatix::adjust_pvalue(method = "fdr")
# Now calculating fold change
log2fc = background_faecium %>%
select(node2,
background_present,
faecium_present) %>%
mutate(n_background_proteins = n_background_proteins,
n_faecium_proteins = n_faecium_proteins,
proportion_faecium= faecium_present/n_faecium_proteins,
proportion_background = background_present/n_background_proteins,
log2fc = log2(proportion_faecium/proportion_background)) %>%
select(node2,log2fc,proportion_in_condition = proportion_faecium,
proportion_in_background = proportion_background)
# all metrics
faecium_metrics = corrected_pvalues %>%
inner_join(log2fc,by = c("cytokine" ="node2")) %>%
mutate(enrichment_score = -log10(p.adj) * log2fc,
condition = "faecium")
write_lines(up_in_faecium,"01_r_processed_data/cytokine_inference/faecium_protein_list")D - Plot of stats
all_metrics = dplyr::bind_rows(infected_metrics,healthy_metrics,faecalis_metrics,faecium_metrics) %>%
mutate(direction_of_change = case_when(log2fc > 0 ~ "increased",
TRUE ~ "decreased"),
)
all_metrics$condition = factor(all_metrics$condition, levels=c("infected","healthy","faecalis","faecium"))
s8d = all_metrics %>%
ggplot(aes(reorder(cytokine, -p.adj),p.adj,color = condition,size = direction_of_change))+
geom_point()+
coord_flip()+
geom_hline(yintercept = 0.05,linetype = "dashed",color = "red")+
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90))+
xlab("Cytokine")+
scale_colour_manual(values = c("#FFA500","#FDE725FF","#440154FF","#2A788EFF"))
s8d
ggsave("03_supplemental/s8_cytokine_inference/s8d.pdf",plot = s8d,width = 5,height = 12,units = "in")E - Network of infected
Figure plotted via cytoscape.
# Since we have
cytokines_in_infected_network = infected_metrics %>%
filter(p < 0.05) %>%
pull(cytokine)
#writing a file so we can look at the network further.
infected_cytokine_network = string_db_results %>%
filter(node2 %in% cytokines_in_infected_network,
`#node1` %in% up_in_infected)
write_csv(infected_cytokine_network,"01_r_processed_data/cytokine_inference/03_infected_protein_cytokine_network_filtered.csv")S9 - Prediction of Clinical Outcome by Nanopore Sequencing
A - Sequencing Data from Strains of faecalis, colored by clinical outcome
## Loading in all the data ###
md = read_csv("00_r_input_data/clinical_metadata.csv")
gene_data = read_csv("00_r_input_data/gene_presence_absence_long.csv")
## Combining data and metadata ##
gene_combined = inner_join(md,gene_data,by = "sample_id")
# Formating ##
faecalis = gene_combined %>%
filter(condition == "faecalis") %>%
filter(!is.na(death_during_admission)) %>%
select(sample_id,death_during_admission,Gene,value) %>%
pivot_wider(names_from = Gene,values_from = value)
## Getting into a format where we can do the clustering ##
faecalis_data = faecalis[,3:28290] %>%
#Remove columns that only contain 0s
select_if(colSums(.) != 0) %>%
as.data.frame()
## Fixing the row names so they match for the heatmap ##
rownames(faecalis_data) = faecalis$sample_id
faecalis_annotation = faecalis[,2] %>%
as.data.frame()
rownames(faecalis_annotation) = faecalis$sample_id
## Generating a dendrogram without a heatmap ##
library(dendextend)
d = dist(faecalis_data,method = "binary")
dend = as.dendrogram(hclust(d, method = "ward.D2"))
# let's add some color: Red = dead, black = okay
colors_to_use = c("#000000","#FF0000")[as.numeric(as.factor(faecalis_annotation$death_during_admission))]
#getting in the same order as dendrogram
colors_to_use = colors_to_use[order.dendrogram(dend)]
# Now we can color the brancehes
dend = branches_color(dend,col = colors_to_use,) %>%
set("branches_lwd", 8)
library(pheatmap)
## Generating a dendogram with a heatmap ##
faecalis_annotation = faecalis_annotation %>%
mutate(death_during_admission = as.character(death_during_admission)) %>%
select(mortality = death_during_admission)
#plotting with a heatmap
p = pheatmap(faecalis_data,
clustering_distance_rows = "binary",
clustering_method = "ward.D2",
annotation_row = faecalis_annotation,
labels_col=rep("",ncol(faecalis_data)))
save_pheatmap_pdf <- function(x, filename, width=12, height=8) {
pdf(filename, width = width, height = height)
grid::grid.newpage()
grid::grid.draw(x$gtable)
dev.off()
}
save_pheatmap_png <- function(x, filename, width=12, height=8,res = 150) {
png(filename, width = width, height = height,res = 150,units = "in")
grid::grid.newpage()
grid::grid.draw(x$gtable)
dev.off()
}
#Have trouble manipulating the pdf in inkscape
save_pheatmap_pdf(p,"03_supplemental/s9_sequencing_based_prediction_of_clinical_outcome/s9a.pdf")## png
## 2save_pheatmap_png(p,"03_supplemental/s9_sequencing_based_prediction_of_clinical_outcome/s9a.png")## png
## 2B - Sequencing Data from Faecium strains, colored by clinical outcome
## Processing the data to only contain the faecium isolates ###
faecium = gene_combined %>%
filter(condition == "faecium") %>%
filter(!is.na(death_during_admission)) %>%
select(sample_id,death_during_admission,Gene,value) %>%
pivot_wider(names_from = Gene,values_from = value) %>%
# Removing weird sample that was found to be E.coli.
filter(sample_id != "S68")
# Processing to get rid of columns we dont need
faecium_data = faecium[,3:28290] %>%
#Remove columns that only contain 0s
select_if(colSums(.) != 0) %>%
as.data.frame()
# Changing the row names so we can plot
rownames(faecium_data) = faecium$sample_id
# Extracting the annotations
faecium_annotation = faecium[,2] %>%
as.data.frame() %>%
mutate(death_during_admission = as.character(death_during_admission)) %>%
select(mortality = death_during_admission)
rownames(faecium_annotation) = faecium$sample_id
p2 = pheatmap(faecium_data,
clustering_distance_rows = "binary",
clustering_method = "ward.D2",
annotation_row = faecium_annotation,
labels_col=rep("",ncol(faecium_data)))
save_pheatmap_pdf(p2,"03_supplemental/s9_sequencing_based_prediction_of_clinical_outcome/s9b.pdf")## png
## 2save_pheatmap_png(p2,"03_supplemental/s9_sequencing_based_prediction_of_clinical_outcome/s9b.png")## png
## 2S10 - Checking Confounding Variables
Faecalis vs Faecium
A- IG abundances between patients with no transplant
md = read_csv("00_r_input_data/clinical_metadata.csv")
protein_sum = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
anot = read_csv("00_r_input_data/all_protein_annotations.csv")
sum_annotated = left_join(protein_sum,anot,by = c("uniprotid" = "Entry")) %>%
left_join(md,by = c("BioReplicate" = "sample_id")) %>%
filter(transplant_type == "none") %>%
#Removing one problematic immunglobulin. Negative values for 3 of the samples messed things up.
#Doesn't change the story, just the plot.
filter(uniprotid != "A0A0B4J2B5")
ig = sum_annotated %>%
filter(grepl(".*Immunoglobulin.*",Protein.names)) %>%
filter(!grepl(".*receptor.*",Protein.names))
stat = ig %>%
rstatix::t_test(Abundance ~ Condition) %>%
rstatix::add_significance()
s10a = ig %>%
ggplot(aes(Condition,Abundance))+
geom_violin()+
geom_boxplot()+
ggpubr::stat_pvalue_manual(stat,
y.position = 22,
step.increase = 0.1)+
theme_classic()
s10a
ggsave("03_supplemental/s10_checking_confounding_variables/s10a.pdf",s10a)B- Between Smoking and Non Smoking
smoking_associated = c("APOC3","IGHV3-7","IGLV1-47","RBP4")
md = read_csv("00_r_input_data/clinical_metadata.csv")
protein_sum = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
anot = read_csv("00_r_input_data/all_protein_annotations.csv")
sum_annotated = left_join(protein_sum,anot,by = c("uniprotid" = "Entry")) %>%
left_join(md,by = c("BioReplicate" = "sample_id")) %>%
filter(Gene %in% smoking_associated,
condition != "healthy",
smoking == FALSE)
plotlist = list()
for(i in smoking_associated){
f = filter(sum_annotated,Gene %in% i)
stat = f %>%
rstatix::t_test(Abundance ~ condition) %>%
rstatix::add_significance() %>%
rstatix::add_xy_position()
p1 = ggplot(f,aes(condition,Abundance)) +
geom_violin()+
geom_point()+
ggtitle(i)+
ggpubr::stat_pvalue_manual(stat)
plotlist[[i]] = p1
}
library(gridExtra)
s10b = do.call("grid.arrange", plotlist)
ggsave("03_supplemental/s10_checking_confounding_variables/s10b.pdf",s10b)####C- Remaining Transplant Type associations
transplant_type_associated = c("APOC3","AZGP1","PCYOX1","RBP4","SERPINC1")
md = read_csv("00_r_input_data/clinical_metadata.csv")
protein_sum = read_csv("01_r_processed_data/quant_msstats_protein_level_data.csv")
anot = read_csv("00_r_input_data/all_protein_annotations.csv")
sum_annotated = left_join(protein_sum,anot,by = c("uniprotid" = "Entry")) %>%
left_join(md,by = c("BioReplicate" = "sample_id")) %>%
filter(Gene %in% transplant_type_associated,
transplant_type == "none")
plotlist = list()
for(i in transplant_type_associated){
f = filter(sum_annotated,Gene %in% i)
stat = f %>%
rstatix::t_test(Abundance ~ condition,) %>%
#rstatix::add_significance() %>%
rstatix::add_xy_position(fun = "median")
p1 = ggplot(f,aes(condition,Abundance)) +
geom_violin()+
geom_point()+
ggtitle(i)+
ggpubr::stat_pvalue_manual(stat,)
plotlist[[i]] = p1
}
library(gridExtra)
s10c = do.call("grid.arrange", plotlist,)
ggsave("03_supplemental/s10_checking_confounding_variables/s10c.pdf",s10c)D- Metabolite Associations
metabolomics_data = read_csv("01_r_processed_data/normalized_metabolomics.csv")
gnps_annotations = read_csv("00_r_input_data/gnps_annotations.csv")
m = inner_join(metabolomics_data,gnps_annotations,by = "row_id") %>%
filter(condition != "Healthy",
transplant_type == "none")
topn = read_csv("01_r_processed_data/top10_faecalis_v_faecium_metabolomics_biomarkers.csv") %>%
pull(row_id)
transplant_type_associated = c("RID_21931","RID_22049","RID_22413","RID_22552","RID_11530","RID_26935")
plotlist = list()
for(i in transplant_type_associated){
f = filter(m,row_id %in% i)
name= unique(f$compound_name_cleaned)
stat = f %>%
rstatix::t_test(norm_value ~ condition) %>%
rstatix::add_significance() %>%
rstatix::add_xy_position()
p1 = ggplot(f,aes(condition,norm_value)) +
geom_violin()+
geom_point()+
ggtitle(paste0(name))+
ggpubr::stat_pvalue_manual(stat)
plotlist[[i]] = p1
}
s10d = do.call("grid.arrange", plotlist,)
ggsave("03_supplemental/s10_checking_confounding_variables/s10d.pdf",s10d)
m